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1.  INTRODUCTION 


The  earth's  gravitational  potential,  V,  at  a  point  with  geocentric  coordinates 
(r,  <t> ' ,  x),  may  be  represented  by  a  spherical  harmonic  expansion  (Moritz,  1980a, 
Sec.  3): 


V-  V(r,  *', 


(1.1) 


where  GM  is  the  gravitational  constant  times  the  earth's  mass,  a  is  the  equatorial 
radius  of  the  best-fitting  earth  ellipsoid;  (C^m,  Snm)  and  Pnm  (function  of 
sin  0',  <p‘  is  geocentric  latitude)  are  respectively  the  fully  normalized  poten¬ 
tial,  or  harmonic,  coefficients  and  associated  Legendre  functions  of  degree  n 
and  order  m.  If  we  subtract  the  potential  coefficients  corresponding  to  the 
earth‘d  normal  gravity  field  from  (Cnm,  Snm),  the  residual  potential  coefficients 
(Cnm>  Snm)  when  ubstituted  in  (1.1)  will  yield  the  anomalous  potential  T. 

As  other  quantities  related  to  the  earth's  gravity  field,  e.g.,  gravity 
anomalies,  geoid  undulations,  deflections  of  vertical,  etc.,  are  simply  related 
to  the  anomalous  potential  T,  a  knowledge  of  the  potential  coefficients  (Cnm,  Snm) 
serves  to  describe  the  earth's  gravity  field.  The  computation  of  gravity  related 
quantities  from  the  potential  coefficients  may  be  termed  as  'harmonic  synthesis'. 
Several  very  efficient  algorithms  are  now  available  for  harmonic  synthesis.  A 
comparison  of  these  algorithms  is  described  by  Tscherning  et  al .  (1983).  The 
inverse  operation  of  computing  potential  coefficients  from  a  global  data  set  of 
gravity  related  quantity  may  be  termed  as  'harmonic  analysis'.  The  algorithms 
for  harmonic  analysis  will  be  described  in  Sections  2  and  3  of  this  study  following 
the  development  by  Colombo  (1981). 

The  resolution  of  gravity  related  quantities  computed  by  harmonic  synthesis 
depends  on  the  degree  n  up  to  which  the  potential  coefficients  are  available. 

A  dramatic  improvement  in  the  resolution  of  geoid  undulations  around  Japan  com¬ 
puted  with  potential  coefficients  up  to  degree  n=180,  when  compared  with  n=36, 
may  be  seen  in  Rapp  (1982,  Figures  1  and  2).  Another  advantage  of  potential 
coefficient  determination  to  a  high  degree,  through  harmonic  analysis,  is  the 
use  of  residual  quantities,  referred  to  this  high  degree  field,  for  evaluating 
integral  formulas,  e.g.  computing  geoid  undulations  from  residual  gravity  anom¬ 
alies.  A  much  smaller  residual  data  'cap'  would  be  adequate  with  the  use  of 
high  degree  field.  For  further  discussion  of  the  role  of  high  degree  field  in 
solving  geodetic  problems,  see  Tscherning  (1983). 

A  high  degree  potential  coefficient  field  GEM10C  to  degree  180  was  developed 
by  Lerch  et  al .  (1981)  using  global  l°x  1°  mean  geoid  heights  H*,  residual  to 
GEM  10  B  (complete  to  degree  36),  by  implementing: 


cos  mx 


sin  mx 
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Rapp  has  developed  several  high  degree  fields,  the  current  one  to  degree  180 
is  described  in  Rapp  (1981).  He  utilized  global  l°x  1°  mean  gravity  anomalies 
which  were  first  adjusted  with  a  potential  coefficients  set  to  degree  36.  The 
procedure  will  be  discussed  more  fully  in  Section  2.3.  In  concept,  the  following 
relation  was  implemented: 


(c*  ) 

'  nm  f  1  r  r 

T  /  4TrY(n-l) 

I  S*  V  ° 

f  nm  } 


Ag 


Pnm(sin  <j)‘ ) 
nm 


!cos  mx 
sin  mx 


da 


(1.3) 


where  Ag  is  the  gravity  anomaly  over  a  block  size  da,  Y  is  an  average  value 
of  gravity  over  the  globe,  denoted  by  a,  and  other  notations  are  as  in  (1.1). 

The  main  problems  in  implementing  (1.2)  or  (1.3),  or  similar  relations,  are 
as  follows.  Firstly,  the  uncertainty  in  the  gravity  related  quantity,  H*,  Ag, 
etc.,  is  not  taken  into  account  in  computing  the  potential  coefficients.  These 
coefficients  remain  the  same  so  long  as  the  global  data  set,  e.g.  of  gravity 
anomalies,  is  not  changed,  whether  the  uncertainty  ascribed  to  the  anomalies  is 
20  mgals,  or  5  mgals,  globally,  or  more  realistically  varies  over  the  globe. 
Secondly,  the  error  estimates  of  the  potential  coefficients  are  not  obtained 
satisfactorily  taking  into  account  both  the  uncertainty  or  'noise'  of  the  gravity 
related  quantity,  and  the  finite  block  size,  or  'sampling',  over  which  this 
quantity  is  given  as  a  mean  value.  Thirdly,  in  the  practical  implementation  of 
(1.2),  (1.3),  etc.,  the  integration  is  replaced  by  a  summation,  and  the  fact 
that  we  are  using  mean  values  instead  of  point  values  requires  a  choice  of  'de¬ 
smoothing'  factors  (see  Rapp  (1981),  p.  4),  depending  on  the  degree  of  potential 
coefficients  being  computed.  If  a  few  desmoothing  factors  are  chosen,  separately 
for  different  degree  bands,  to  span  the  entire  spectrum,  a  sharp  discontinuity 
may  result  at  the  boundary  of  degree  bands.  In  fact,  the  current  potential  co¬ 
efficients  field  by  Rapp  (1981)  was  developed  to  degree  300,  but  has  been  usually 
employed  only  up  to  degree  180  because  of  concern  over  such  a  discontinuity. 

Colombo  (1981)  developed  algorithms  for  'optimal'  estimation  of  potential 
coefficients  which  are  free  of  the  above  problems.  His  tests  were  carried  out 
with  global  set  of  5°x  5°  anomalies.  He  estimated  CPU  time  of  about  two  hours, 
on  an  Amdahl  470  V/6-II  computer,  for  computing  coefficients  to  degree  180  with 
a  global  set  of  l°x  1°  anomalies.  With  a  faster  470  V/8  computer,  and  by  slight 
modifications  in  the  subroutines,  we  can  now  compute  coefficients  to  degree  250 
with  l°x  1°  anomalies  in  about  60  minutes  CPU  time.  This  procedure  will  be 
described  in  Section  3. 

The  initial  potential  coefficients  set  used  in  this  study  is  the  one  developed 
by  Rapp  (1981)  to  degree  180.  We  used  this  set  to  compute  a  global  l°x  1°  anomaly 
data,  which  then  served  as  a  test  data  to  generate  different  potential  coefficient 
sets  by  optimal  estimation  procedure,  assuming  different  error  estimates  for  the 
anomaly  data.  Rapp  (ibid.)  had  also  generated  an  adjusted  global  l°x  1°  anomaly 
data  set,  adjusted  with  a  potential  coefficient  set  to  degree  36.  This  was  the 
second  global  anomaly  data  set  used  in  the  present  study.  The  different  poten¬ 
tial  coefficient  sets,  and  the  global  anomlay  data  sets,  are  described  in  Section  4. 
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The  differences  between  the  various  potential  coefficient  sets  are  examined 
i  Section  5.  This  allows  us  to  estimate  the  improvement  by  the  optimal  estimation 
rocedure  over  the  current  (Rapp,  1981)  harmonic  analysis  procedure  in  developing 
high  degree  field.  There  are  presently  several  areas  in  the  world,  which  have 
aophysically  predicted  l°x  1°  anomalies.  We  therefore  also  examine  in  Section 
.3,  the  change  in  potential  coefficients  set  if  anomalies  in  some  areas  of  the 
orld  are  set  to  zero. 

The  potential  coefficients  developed  to  degree  and  order  250,  and  their 
rror  estimates,  by  optimal  estimation  procedure  are  discussed  in  Section  6. 
here  is  a  very  substantial  improvement  in  these  error  estimates  when  compared 
ith  the  error  estimates  of  the  current  set  of  potential  coefficients  (Rapp, 

981).  The  error  estimates  of  geoid  undulations,  and  gravity  anomalies,  computed 
rom  the  potential  coefficients  developed  by  optimal  estimation,  are  also 
xamined  by  degree,  and  cumulatively.  Section  7  summarizes  the  findings  of  the 
resent  study.  The  choice  of  the  highest  degree  250,  to  which  the  coefficients 
'ill  be  estimated  from  l°x  1°  anomalies,  was  somewhat  arbitrary.  It  was  chosen 
rom  the  dual  consideration  of  keeping  the  computational  effort  manageable,  and 
et  not  lose  the  high  degree  information  available  in  l°x  1°  anomalies.  The  re- 
ormulation  of  the  computational  algorithm  in  several  steps  is  given  later  in 
able  3.4.  We  will  also  find  in  Section  6,  that  at  degree  250,  the  undulation 
lagnitude  per  degree  is  about  1  cm,  and  the  anomaly  magnitude  per  degree  is  about 
1.5  mgals. 


ut  on  tape  after  rearranging  as  180  ( =Ne (2.22) )  numbers,  (NN+1)  times,  for  90 
atitude  bands.  The  requirement  is  to  form  them  into  (NN+l)(R(m),  m=0  to  NN) 
atrices  of  180  x  180  numbers.  The  expansion  of  180  x  90  numbers  into  180x180 
atrix,  by  exploiting  the  persymmetric  nature  of  R(m)  matrices,  can  be  handled 
y  proper  indexing.  However,  reading  180  numbers  as  one  logical  record  still 
equires  90  x  (NN+1)  records  to  be  read  to  access  the  whole  set  to  pick  out  elements 
or  one  (R(m)  matrix;  and  90 x  (NN+1)2  records  to  be  read  to  form  (NN+1)  R(m) 
atrices,  and  additional  TAPE  10  for  writing  out  the  reordered  R(m)  matrices, 
sing  the  largest  block  size  for  reading  and  writing,  the  TAPE  10  for  different 
N  are  given  in  Table  3.2. 


able  3.2  TAPE  10  for  Reordering  Normals  for  Estimating  Coefficients  to  Degree  NN. 


NN 

TAPE  10 

Read 

Write 

Total 

180 

134,  121 

724 

134,  845 

250 

257,  777 

1,004 

258,  781 

300 

370,  832 

1,204 

372,  036 

To  avoid  excessive  TAPE  10  costs,  the  coding  was  modified  to  put  the  sequen- 
:ial  set  of  generated  normals,  after  grouping  180  x  (NN+1 )  numbers  for  k  latitude 
>ands  as  one  record,  on  random  disk  storage.  This  not  only  removed  the  need  to 
'ead  the  entire  sequential  set  to  form  one  R(m)  matrix,  but  also  allowed  k  R(m) 
latrices  to  be  assembled  at  a  time,  as  any  one  record  was  read  in  from  random 
itorage.  However,  this  required  that  an  array  of  kx  (16,290)  x8  bytes  should  be 
ivailable,  where  16,290  (  =  180  x  181/2)  is  the  number  of  double  precision  elements 
in  vector  form  of  a  180  x  180  R(m)  matrix. 

With  2048  K  bytes  of  virtual  storage,  k  was  chosen  as  13.  A  proportionally 
larger  value  of  k  could  be  chosen  if  a  larger  region  size  is  available.  This 
;tep  of  reordering  the  normals  may  be  termed  as  step  3A. 

1.32  Reordering  of  the  Right  Hand  Side  Vectors 

After  the  reordering  of  R(m)  matrices,  m=0  to  NN,  we  need  to  read  in  the 
'ight  hand  side  vector  in  the  following  manner,  for  the  solution  of  2Lnm  in 
(2.30).  For  j£r,m,  180  elements  for  latitude  bands  i=0  to  179,  arranged  by  order 
i=0  to  NN,  and  degree  n=m  to  NN;  where  NN  is  the  highest  degree  and  order  up 
;o  which  coefficients  are  to  be  estimated,  usually  N<NN<Nmax.  From  the  45,451 
[i  values  per  latitude  band,  n=0  to  Nmax,  m=0  to  n,  available  in  step  1C,  the 
/aiues  corresponding  to  m=0  to  NN,  n=m  to  NN,  could  be  picked  by  suitable  indexing 
in  each  latitude  band,  and  multiplied  by  the  factor  in  (3.4)  to  obtain  corres- 
jonding  kjL.  This  would  result  in  the  elements  of  k^  being  assembled  as 
(NN+1 )(NN+2)/2  elements,  m=0  to  NN,  n=m  to  NN,  for  each  latitude  band  i=0  to 
179.  If  we  consider  it  as  a  general  matrix  stored  in  a  vector  form  row  by  row, 
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m=0  to  NN,  for  each  row  p.  The  procedure  is  repeated  by  reading  in  1^  for 
i  =  l,  2,  ....  89  and  implementing  p=i  to  179-i  for  each  i,  thereby  fully  exploiting 
the  equatorial  symmetry  of  the  data  grid,  which  makes  R(m)  matrices  persymmetrical , 
i.e.,  symmetrical  wi  h  respect  to  both  the  main  diagonal  and  the  main  antidiagonal. 

For  a  given  NN,  the  CPU  time  for  generating  r1^  varies  almost  linearly  with 
the  latitude  band  i,  decreasing  from  high  to  low  latitudes,  except  for  slight 
reduction  in  the  polar  region  for  high  NN.  A  modification  was  made  to  generate 
only  a  portion  of  the  normals,  i.e.  HP  for  any  one  or  more  i  (for  all  p=l  to 
179-i)  for  estimating  the  time  to  generate  the  entire  normals  for  a  given  NN. 

Some  CPU  times  on  an  Amdahl  470  V/8  computer  are  gr'en  in  Table  3.1  for  various 
NN  for  the  case  of  l°x  1°  anomalies. 


Table  3.1  CPU  Time  for  Generating  Normals  for  Estimating  Coefficients  to  Degree 


Latitude  Band 


CPU  Time  in  Seconds 


i 

NN=12 

NN=60 

NN=250 

0 

5.3 

13.4 

25.7 

29 

4.1 

8.1 

31.8 

59 

2.9 

5.0 

10.0 

89 

1.7 

1.8 

2.1 

Total  Time 

=  5.3  min. 

=  10  min. 

=  26  min. 

i=0  to  89 

The  generation  of  normals  in  this  step  2  depends  only  on  the  data  grid, 
e.g.  Ae=AA=l°.  As  the  error  estimates  of  anomalies  are  introduced  in  later  steps, 
the  normals  are  not  required  to  be  regenerated  for  testing  the  effect  of  different 
anomaly  error  estimates. 

3.3  Optimal  Quadrature  Weights  and  Error  Estimates 

The  elements  of  R(m)  matrices  generated  in  step  2  in  Section  3.2  first  need 
to  be  arranged  by  the  order  m,  m=0  to  NN,  in  a  vector  form  for  the . 180  x  180  symmetri 
matrix,  for  all  latitude  bands  i  and  p.  Similarly,  the  elements  llm  in  step  1C 
need  to  be  arragned  by  the  order  m  to  assemble  the  vector  J<nm  in  tne  form  needed 
for  implementing  (2.30)  to  solve  for  Xnm>  for  all  n  for  a  given  m.  It  is  only 
at  this  stage  that  average  variance  o?  of  anomalies  in  (2.32)  in  different  latitude 
bands  has  to  be  considered. 


3.31  Reordering  of  the  Normals 

r1^  are  generated  in  step  2  as  (NN+1)  numbers,  for  180  latitude  bands  p, 
for  90  latitude  bands  i  from  pole  to  equator.  The  generated  numbers  are  written 


even 


(3.3) 


PI 


i 

nm’ 


(i  =  0  to  89) 


179-i ) ;  n  -  m  = 

odd 


Finally,  as  the  PInm  are  needed  in  the  form  of  Ipm 
further  computations,  a  third  tape  was  written  for  O 
pole,  i  =  0  to  179: 


as  in  (2.24)  for  all 
from  north  pole  to  south 


ri  =  /  cn  .  _L  .  PI i  = 
nm  \2n+l  a  .  nm 


(3.4) 


These  three  steps  of  computing  Plpm*  1  *  0  to  89;  extending  PIpm  to  i  =0  to  179; 
and  modifying  these  to  I^m,  i  =  0  to  179  are  preliminary  steps  to  the  generation  of 
normals.  These  three  steps  may  be  termed  as  1A,  IB  and  1C  for  convenience  of  refer¬ 
ence.  The  anomaly  degree  variance  model  will  be  specified  in  Section  4. 

3.2  Generation  of  Normals 

The  Fourier  transforms,  a {f,  m  =  0  to  n,  of  the  covariances  of  anomaly  blocks 
in  two  latitude  bands  i  and  p  are  obtained  by  (Colombo,  1981,  p.  44,  (2.63)): 


a 


ip 

m 


Nmax 


-  [  l 

n=m 


I 


i 

n,2N-m 


I 


P 

n,2N-m 


]  * 


AX2  for  m=0 

2 

^2- (1-cos  m  ax)  for  m^O 


(3.5) 


where  2N-m  }  Nmax.  The  second  term  in  the  square  brackets  take  the  aliasing  of 
Fourier  transform  into  account,  with  increasing  m.  Nmax  was  taken  as  300  based 
on  experiments  by  Colombo  (1981,  Sec.  4.2). 

The  element  r1^,  corresponding  to  latitude  bands  i  and  p,  of  the  matrix 

R(m),  m  =  0  to  NN,  where  NN  is  the  highest  degree  up  to  which  coefficients  will 

be  estimated  (N  <  NN  <  N  „  ),  are  obtained  by: 

max 


| 2N,  if  m=0 
In  ,  if  m?0 


(3.6) 


The  formation  of  R(m)  matrix  may  be  termed  as  the  generation  of  normals 
in  the  frequency  domain  corresponding  to  the  formation  of  the  covariance  matrix 
CZ2  of  the  64,800  anomaly  blocks.  This  step  of  the  algorithm  takes  a  large 
CPU  time.  However,  the  CPU  time  is  dependent  on  the  highest  degree  up  to  which 
the  coefficients  are  to  be  estimated,  as  the  number  of  180  x  180  R(m)  matrices 
to  be  formed  is  NN+1. 

The  implementation  of  (3.5)  starts  by  reading  in  the  45,451  values  of  Inm 
in  (3.4)  in  step  1C  for  i=0,  and  multiplying  with  the  corresponding  values  of 
iRm  1n  (3.5)  for  p=0  to  179  yielding  NN+1  numbers  a°P  (and  r°P  from  (3.6)), 
nm  mm 


3.  ALGORITHM  FOR  OPTIMAL  ESTIMATION  OF  POTENTIAL  COEFFICIENTS 


The  algorithm  outlined  in  Section  2.2  was  coded  by  Colombo  (1981,  Appx.  B4) 
in  subroutines  NORMAL  and  ANALYS,  and  were  tested  by  him  up  to  the  smallest 
anomaly  block  size  of  5°x5°,  comprising  of  2592  anomalies  globally.  In  applying 
these  subroutines  to  the  case  of  64,800  l°x  1°  anomalies,  it  appeared  advisable 
both  for  testing,  and  for  actual  computations,  to  separate  out  portions  of  NORMAL 
which  could  run  sequentially,  without  repeating  earlier  stages  and  thus  allowing 
greater  flexibility  in  considering  different  anomaly  uncertainties,  and  later 
different  global  anomaly  data  sets.  Also,  two  major  modifications  were  required 
because  of  very  large  TAPE/DISK  10' s,  or  very  large  array  size,  needed  with 
the  existing  coding  for  the  case  of  l°x  1°  anomalies.  The  separate  portions  of 
implementing  the  algorithm  and  the  modifications,  are  described  below.  The 
numerical  values,  and  ranges,  will  refer  to  the  case  of  l°x  1°  anomalies. 

3.1  Integration  of  Associated  Legendre  Functions 

The  integrals  of  associated  Legendre  functions,  PI^  , 


PI  =  f  1  P  (cos  e)  sinedo 
nm  J  nm 

ei 

are  needed  repeatedly  for  implementing  (2.20)  and  (2.21),  e.g.  in  the  formation 
of  the  elements  of  R(m)  matrices  and  the  array  for  computing  the  optimal 
quadrature  weights  Xj,m  in  (2.30).  The  recursive  relations  for  efficiently  computing 
PI^  were  developed  by  Paul  (1978).  A  subroutine  supplied  by  him  (private  com¬ 
munication  to  R.H.  Rapp)  was  used  to  generate  £1  on  tape  for: 


0°  <  0.  <  89°,  A0  =  1°,  n  =  0  to  N„  :  m  =  0  to  n;  N  =  300 
1  max  max 


(3.2) 


These  are  nominal  values  of  Ae,  and  0  as  the  complement  of  geodetic  latitude 
ip.  The  actual  integration  limits  in  (3.1)  were  obtained  by  0  =  90  -  <p ' ,  by  con¬ 
verting  the  nominal  geodetic  latitudes  <j>  to  the  corresponding  geocentric  latitudes 
<f>‘.  To  avoid  the  numerical  instability  in  computing  the  sectorial  PljL,  the 
direct  series  expansions  were  used  for  latitudes  up  to  50°  from  the  pole  for  PT 
instead  of  forward  sectorial  recurrence  relations.  This  tape  was  already  avail¬ 
able  for  this  study,  but  for  any  future  studies  with  n  exceeding  300,  Gleason 
(1983)  has  suggested  modifications  to  Paul's  routines  for  employing  backward  and 
forward  sectorial  recursions  depending  on  the  condition  number  associated  with 
the  recursion  (Gerstl,  1980). 

The  order  of  _PI  on  tape  was  according  to  latitude  bands  from  the  north  pole 
to  the  equator;  and  on  each  latitude  band  (Nmax+l)(Nmax+2)/2=45,451  P_P  values 
were  arranged  by  degree  and  order  as  in  (3.2).  As  the  PIJL  are  required  for  all 
latitude  bands  from  the  north  pole  to  the  south  pole  for  tne  elements  in  the 
right. hand  side  of  (2.30),  a  new  tape  was  written  by  exploiting  the  symmetry 
of  P I Im  with  respect  to  the  poles,  and  change  of  sign  in  the  two  hemispheres 
when  (n-m)  is  odd: 


We  note  from  (2.40)  that  under  the  assumption  of  a  constant  global  anomaly 

error,  the  propagated  error  estimate  of  all  coefficients  of  different  order  are 

the  same  in  each  degree.  This  also  appears  to  be  the  case  in  (2.41)  for  sampling 

error,  but  this  is  an  artifact  due  to  the  fitting  of  the  quartic  expression 

in  (2.41)  to  estimation  errors  0en  per  degree  in  (2.33)  instead  of  estimation 

errors  for  each  pair  of  coefficients  C“  as  discussed  after  (2.32). 

nm 

Based  on  root  mean  square  errors  of  56,751  merged  anomalies 
as  ±10  mgals,  and  ±30  mgals  for  the  balance  8049  anomalies,  the  constant  global 
anomaly  error  m  in  (2.36)  could  be  taken  as  ±15  mgals.  A  larger  value  of  ±20 
mgals  was  chosen  to  allow  for  possible  anomaly  errors  in  the  several  largely 
unsurveyed  areas  (Rapp,  1981,  p.  15).  This  caused  a  slightly  pessimistic  coef¬ 
ficient  error  estimate  in  (2.40). 


where  a  is  the  standard  error  of  a  l°x  1°  anomaly,  assumed  to  be  the  same  for  all 
64,800  anomaly  blocks;  and  pn  are  the  de-smoothing  factors  (2.7)  in  HARMIN. 

are  the  partials  in  (2.34)  and  P^  is  the  weight  matrix  of  anomalies, 
assumed  to  be  diagonal. 


[P  ]  =  cos  <f>7o2  =  sin  e/a2 


(2.37) 


The  corrections  V_  to  the  anomaly  data  set  is  obtained  by: 


V  =  P~ 1 Bt  P  V 
—l  l  —i  x— x 


(2.38) 


The  simplification  in  the  adjustment  model  results  from  B^P^bJ  matrix  being 
diagonal  (so  that  its  inversion  becomes  trivially  simple),  which  depends  on  the 
assumption  (2.37)  with  the  standard  error  a  of  all  anomalies  being  the  same  globally. 

Elements  of  the  covariance  matrix  Enx  of  coefficient  errors  due  to  propagation 
of  anomaly  errors  are  given  (Rapp,  1969,  p.  112,  (12))  by: 


[E 


nx 


]  - 


A 

1+[PX]A 


(2.39) 


E  was  considered  as  a  diagonal  matrix.  The  standard  errors  of  'adjusted  SET  1  ‘ 
coefficients  to  degree  36  (some  additional  coefficients  to  degree  48)  were  retained, 
while  for  higher  degree  and  order  coefficients  Px  was  null  matrix  and  pn  in  (2.7) 
was  uniformly  taken  as  1/4-n .  Then  the  diagonal  elements  of  Enx  for  coefficients, 
other  than  the  'adjusted  SET1',  were  computed  (Rapp,  1981,  p.  29,  (30))  by  using 
(2.36)  as: 


[E 


nx 


]  - 


a2  A<j>  AA 

4ir[Y(n-l)]2 


(2.40) 


The  diagonal  covariance  matrix  Ens  of  coefficient  errors  due  to  sampling 
finite  size  of  l°x  1°  anomaly  blocks  was  based  on  empirical  relation  derived  from 
a  Montecarlo  approach  (Colombo,  1981,  p.  78,  (3.10): 


[Ens]  =  ^l00  { ( ("16, 19570  ($)  +  30.34506)  (£)  +  40.29588)  (jj)2}  ]  2  (2.41) 


The  Montecarlo  experiments  are  described  in  Colombo  (1981,  Sec.  3.1)  and  the 
sampling  error  computations  were  performed  as  described  in  Section  2.2  of  this 
study  utilizing  noise-free  data. 


2.3  Rapp's  December  1981  Potential  Coefficients 

Rapp's  current  set  of  coefficients  was  developed  (Rapp,  1981)  to  degree 
and  order  300,  but  have  been  generally  employed  to  degree  and  order  180.  The 
procedure  for  computing  these  coefficients  is  reviewed  here.  In  concept,  the 
coefficients  were  computed  as  in  Section  2.1,  but  the  error  estimates  included 
sampling  error  based  on  Section  2.2. 

First,  a  low  degree  potential  coefficients  set  (complete  to  degree  36  with 
some  additional  coefficients  to  degree  48)  was  assembled  by  a  weighted  merging 
of  most  current  sets,  including  resonance  coefficients  of  orders  12  to  15  and 
30.  These  merged  coefficients  were  called  SET1,  and  were  also  assigned  error 
estimates.  Similarly,  a  global  set  of  64,800  l°x  1°  mean  gravity  anomalies  was 
assembled  by  merging  most  current  terrestrial  anomaly  data  with  anomalies  deter¬ 
mined  from  adjusted  Seasat  altimeter  data.  The  number  of  merged  anomalies  was 
56,751  with  root  mean  square  of  standard  errors  of  anomalies  being  ±10  mgals. 

The  remaining  8049  anomalies,  to  complete  the  global  set,  were  assigned  values 
implied  by  SET1  coefficients,  and  assigned  standard  error  of  30  mgals. 

Next,  a  simplified  adjustment  was  made,  whose  procedure  is  outlined  below, 
of  SET1  coefficients  and  64,800  l°x  1°  anomalies.  The  result  of  the  adjustment 
may  be  called  'adjusted  SET1 '  to  degree  36  (some  additional  coefficients  to 
degree  48),  and  64,800  'adjusted  anomalies'.  The  'adjusted  anomalies'  were  used 
to  compute  potential  coefficients  to  degree  300  using  subroutine  HARMIN,  which 
implements  (2.5)  with  de-smoothing  factors  (2.7).  The  'adjusted  anomalies'  and 
the  coefficients  (using  HARMIN)  up  to  degree  180  will  be  termed  test  data  set 
B  in  Section  4.2.  Because  of  the  simplified  adjustment  procedure,  the  low  degree 
coefficients  in  data  set  B  did  not  fully  agree  with  the  'adjusted  SET1'  coef¬ 
ficients,  though  the  difference  was  small  as  discussed  later  in  Section  5.  A 
merged  set  of  coefficients  was  then  formed  using  'adjusted  SET1’  where  available 
(complete  to  degree  36,  with  some  additional  coefficients  to  degree  48),  and 
the  remaining  high  degree  coefficients  as  developed  from  'adjusted  anomalies' 
using  HARMIN.  This  merged  set  up  to  degree  180  are  Rapp's  current  (December 
1981)  potential  coefficients,  and  will  be  part  of  test  data  set  A  in  Section  4.2. 

The  simplified  adjustment  of  combining  SET1  with  global  anomaly  data  is 
based  on  the  model : 


F  =  HjJo  JLJ)  =  0  (2.34) 

where  Lj  are  the  adjusted  SET1  coefficients,  while  L®  are  the  adjusted  anomalies. 
Denoting  one  element  as  [•],  an  element  [Vx]  of  the  correction  vector  V^to 
SET1  coefficients  is  given  by: 


[Vx]  =  -[w]/(l  +  A[Px]) 


(2.35) 


where  w  is  the  (Disclosure  vector  between  the  SET1  coefficients  and  the  computed 
coefficients  from  the  anomalies  using  HARMIN,  Px  is  the  weight  matrix  (assumed 
to  be  diagonal)  for  the  SET1  coefficients,  and  A  is  one  of  the  elements  of  the 
diagonal  matrix  j5  P-1  Bj,  given  by  (Rapp,  1981,  p.  6,  (23)): 

it  it  it 
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(2.31)  allows  harmonic  analysis  through  a  quadrature  formula,  and  may 
be  called  the  optimal  quadrature  weights.  The  optimal  quadrature  weights  are 
derived  on  the  basis  of  minimizing  the  sum  of  both  error  quadratic  terms  due  to 
finite  size  of  anomaly  blocks,  and  the  anomaly  uncertainties.  And  (2.31)  retains 
the  efficiency  in  computation  of  the  coefficients  inherent  in  the  structure  of 
(Czz+D)  discussed  before  (2.24).  This  however  requires  (Colombo,  1981,  Sec.  2.9) 
a  slight  modification  in  the  structure  of  noise  matrix  D.  The  modification  is 
that  the  variances  for  all  anomalies  in  a  latitude  band  be  equated  to  the  average 
value  for  that  band,  i.e. 


2N-1 


(2.32) 


This  modification  causes  a  change  in  estimates  for  propagated  error  (en  in 
(2.15))  for  non-zonal  coefficients. _  But  the  sum  of  variance  of  a  pair  of  non- 
zonal  coefficients  C,fm,  i.e.  (Cnm,  Snm),  is  unaltered.  The  estimated  variance 
of  each  of  the  coefficients,  in  a  pair,  is  then  arbitrarily  made  equal.  This 
is  a  minor  modification  to  retain  the  computational  efficiency  of  not  having 
to  directly  invert  the  (Czz+D)  matrix  of  a  large  size. 

There  is  no  change  in  error  estimates  of  the  zonal  coefficients.  There  is 
also  no  change  in  the  error  variance  per  degree,  o|n,  or  the  error  variance  per 
average  coefficient  per  degree,  o|n: 
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(2.33) 


The  implementation  of  (2.23)  to  (2.29),  leading  to  the  computation  of  optimal 
quadrature  weights  in  (2.30)  is  done  by  subroutine  NORMAL.  The  error  estimates 
for  the  coefficients,  obtained  in  concept  through  (2.21),  are  also  computed  by 
NORMAL,  including  (2.33),  under  the  condition  (2.32).  The  optimal  quadrature 
weights,  x^m*  computed  by  NORMAL,  are  used  in  subroutine  ANALYS  with  the  global 
anomaly  data  for  the  harmonic  analysis  of  coefficients  in  (2.31).  The  coding 
of  NORMAL  and  ANALYS  is  documented  in  Colombo  (1981,  Appx.  B4).  Some  changes 
made  in  these  subroutines  will  be  discussed  in  Section  3. 

One  may  notice  the  similarity  of  (2.5)  and  (2.6),  the  quadrature  formula 
with  de-smoothing  factors  un,  and  (2.31)  read  with  (2.28),  (2.26),  (2.30),  i.e. 
the  quadrature  formula  with  optimal  quadrature  weights  x^.  As  already  mentioned, 
x^m  are  dependent  on  the  anomaly  uncertainties  through  (2.30),  while  pn  computed 
through  (2.7)  and  (2.8)  do  not  take  into  account  the  anomaly  uncertainties. 

Hence,  the  coefficients  computed  in  Section  2.1  will  remain  the  same  whether 
the  global  anomaly  data  has  standard  error,  as  an  example,  of  5  mgals,  or  20  mgals, 
or  more  realistically  varies  over  the  globe.  However,  the  estimated  coefficients 
in  Section  2.2  would  real istical ly  reflect  the  uncertainties  in  the  estimate  of 
anomalies,  subject  to  the  averaging  of  variance  in  a  latitude  band  in  (2.32). 
Section  2.2  coefficients  are  also  optimal  in  the  sense  of  minimizing  Ey  in  (2.16), 
and  satisfactorily  accounting  for  both  Es  and  En,  and  also  do  not  show  any  dis¬ 
continuities  which  are  caused  in  Section  2.1  coefficients  due  to  (2.7). 


-^firna ,  z 


2n+l 


,  6,+A9_  x,+ax  T  T 

[•••  7—  I  1  P__(cos  0)  sine  de  J  J  c“  dx  ...]  , 

+1  A-j  J  nm  J  — m 


(2.25) 


1=0  to  (N0-l) 


•  X  .+AX  j  -r 

=  [...  k^m  J  J  c“T  dx  ...]' 
nm  1  — m 

Aj 


where 


im 


[... 


cos 

sin 


m  j  ax  . . .]  ,  j=0  to  ( Nx~ 1 ) 


(2.26) 


(2.27) 


X  -+AX 
and  /  J 


where  A(m) 


cos] 

1 

"| 

A(m) 

* 

sin' 

|  mx  dx  = 

-B(m) 

cos  mj  ax  + 


B(m) 

A(m) 


AX 

(sinm  AX)/m 


,B(m)= 


(cos  m  AX-l)/m 


sin  mj  ax 


if  m=0 
if  m?0 


(2.28) 


(2.29) 


If  we  collect  kpm  for  all  latitude  bands  i=0  to  (N0-l)  to  define  Jcpm  and 
similarly  collect  x,*m  for  all  latitude  bands  to  define  xnm,  See  below,  i.e.. 


k._  =  [k®  ...  kL  ...  kJJ®"1]7 

— nm  nm  nm  nm 


=  r  0  i  "e-V 

^■nm  *-xnm  ' '  *  xnm  *  ‘  *  xnm  -* 


xhere  ^  *  (Rfml+W)-1  k„m  , 


(2.30) 


W  being  the  diagonal  matrix  consisting  of  average  anomaly  variance  in  each  latitude 
band  (see  (2.32)  below);  then  Colombo  (1981,  p.  43,  (2.61))  shows  that  elements  of 
vector  c  in  (2.20)  may  be  computed  by: 

,  N-l  ,  2N-1 

Pa  -  1  y  v1  y 

^nm  y(n-l)  nm 

where  we  have  substituted  N  for  N0  and  2N  for  in  view  of  (2.4)  and  (2.22). 


_% 


V 


A(m) 

B(m) 

" 

■■•S'? 

cos  mj  ax  + 

sin  mj  ax 

Agi0-  (2.31) 

.  ■  ■  j 

_ 

-B(m) 

A(m) 

J 

I 

•v  -.1 
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The  formation  and  inversion  of  (Czz  +  D)  for  64,800  l°x  1°  anomalies  in  (2.20) 
and  (2.21)  would  be  a  formidable  task.  But  Colombo  (1979)  has  shown  that  under 
general  conditions  fulfilled  by  equi-angular  grids,  each  covariance  matrix  between 
l°x  1°  anomaly  blocks  in  two  latitude  bands  has  Toeplitz  circulant  structure 
(Colombo,  1981,  Sec.  2.10).  The  main  requirement  of  the  same  constant  value 
of  AA  in  each  of  the  N  latitude  bands  in  (2.4): 


aa  =  Aj+1  -  A j ;  j  =  0  to  359  for  all  e ^ ;  i  =  0  to  179 

A0  =  ei+1  -  8.;  A9  =  AA  =  1°;  N0  =  =  180;  Nx  =  ~  =  360  (2.22) 


is  naturally  fulfilled  for  equi-angular  grid  (the  numbers  correspond  to  l°x  1° 
grid). 

The  structure  of  (Czz  +  0)  makes  it  possible  that  its  inversion  can  be  equi¬ 
valently  substituted  (Colombo,  1981,  p.  42,  (2.54))  by  the  much  simpler  inversion 
of  N0  x  N0  matrices  R(m),  m=0  to  if  we  wish  to  solve  only  for  coefficients 
up  to  the  Nyquist  frequency  given  by  N0  =  However,  if  we  were  to  solve 

for  coefficients  up  to  degree  NN  >  Ne,  it  will  require  the  formation  of: 


N  x  N  matrices  R(m),  m=0  to  NN 

U  U 


(2.23) 


However,  their  inversion  is  still  a  relatively  simple  problem  as  their  size  is 
Ne  x  N0.  Also,  the  elements  of  R(m),  which  are  related  to  the  discrete  Fourier 
transform  of  covariances  of  anomaly  blocks  in  a  pair  of  latitude  bands,  can  be 
obtained  directly  without  actually  computing  the  anomaly  covariances.  .This  re¬ 
quires  repeated  manipulation  (Colombo,  1981,  Sec.  2.11)  of  integrals  I1’  of 
associated  Legendre  functions  in  the  form: 


nt 


-4 


2n+l 


i  VAe-  i 

•  Pnt(cos  9)  sine  de  =  knt 

i  ef 


(2.24) 


where  subscript  t  is  related  to  the  order  m,  k^t  is  defined  in  (2.25)  below, 
other  notations  are  as  in  (2.1),  (2.4),  (2-6) ;  and  cn  are  the  anomaly  degree 
variances.  The  computation  of  integrals  in  (2.24)  will  be  discussed  in  Section 
3.1;  the  computation  of  the  elements  of  matrices  R(m)  will  be  discussed  in 
Section  3.2,  and  the  anomaly  degree  variance  model  will  be  specified  in  Section  4. 

The  components  c^ma,  z  the  cross-covariance  matrix  Ccz  in  (2.20),  (2.21) 
are  computed,  for  degree’ n— and  order  m,  by: 
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where  e^  represents  the  propagated  noise,  while  may  be  Interpreted  as  the 
sampling  error,  would  tend  to  zero  If  perfect  data  (without  any  error  or 

noise)  was  employed  in  (2.13),  and  the  number  of  data  points  tended  to  infinity; 
or  in  other  words,  the  mean  anomaly  block  size  tend  to  zero.  We  assume  that  the 
'errors  es  and  e^,  arise  due  to  completely  independent  causes. 

The  total  covariance  matrix  Ej  of  the  potential  coefficient  errors,  may 
then  be  written  as  the  sum  of  covariance  matrices  of  errors  Es  due  to  sampling 
and  En  due  to  data  noise: 


Es  +  E„  -  H  (e^J  ♦  M  (e^} 


=  M  {(c.  -  FzJ(c  -  Fz)T}  +  M  {Fnji^F1}  (2.16) 

where  M  {  }  is  the  averaging  operator. 

(2.16)  is  easily  simplified  to: 

ET  =  C  -  2CczFT  +  F(Czz  +  D)  FT  (2,17) 

where 


C  *  M  {ccT},  Cc2  -  M  {czT},  Czz  =  M  {Z£T}  (2.18) 

D  =  M  {n  nT} 


are  the  covariance  matrix  of  the  coefficients,  cross  covariance  matrix  of  coef¬ 
ficients  and  mean  anomalies,  covariance  matrix  of  mean  anomalies,  and  covariance 
matrix  of  mean  anomaly  errors,  respectively. 

An  optimal  value  of  estimator  operator  F  in  (2.13)  may  now  be  obtained, 
which  will  minimize  the  sum  of  square  of  total  errors  for  all  coefficients,  i.e. 
the  trace  of  Ej,  by  setting  3tr[Ej]/aF  to  zero.  This  leads  (Colombo,  1981, 
Sec.  2.8)  to  the  well-known  relations  in  least  squares  collocation: 


F  =  C  (C  +  D)_1 
czv  zz  ' 


(2.19) 


and  by  inserting  (2.19)  in  (2.13)  and  (2.17),  we  get: 


For  the  case  of  harmonic  synthesis,  we  may  write  for  the  point  and  mean 
anomalies: 


However,  as  the  potential  coefficients  are  available  only  up  to  degree 

n  <_  Nmax,  in  the  practical  implementation  of  (2.10)  the  estimate  ^7. .  is  affected 

by  the  truncation  error  at  N  „  (and  similarly  for  (2.9)): 

max 


Nmax 

l 

n=2 


(n-1) 


4 

ij 


(2.11) 


Subroutine  HARMIN  implements  the  harmonic  analysis  (2.5)  using  the  de-smoothing 
factors  in  (2.7).  Subroutine  SSYNTH  implements  the  harmonic  synthesis  (2.11). 

The  coding  of  these  subroutines  is  documented  in  Colombo  (1981,  Appx.  B2,  B3). 

2.2  Quadrature  Formulas  with  Optimal  Weights 

If  we  consider  the  global  mean  gravity  anomaly  vector  Ag  to  comprise  two 
parts,  the  signal  z  and  the  noise  n,  i.e.  , 


Aq  =  z  +  n 


(2.12) 


and  denoting  the  potential  coefficients  vector  [C,fm],  m=0  to  n,  n=0  to  Nmax  by 
c_,  which  is  estimated  from  the  anomaly  data  through  a  system  of  linear  equations 
denoted  by  the  matrix  operator  F,  we  may  write: 


c  =  £(z  +  n) 


(2.13) 


An  optimal  value  of  the  estimator  operator  F  will  be  derived  later  in  (2.19). 

The  estimation  error  vector  e  may  then  also  be  broken  into  two  parts  e 

and  e  :  s 

— n 

£  =  £-  £  =  £-  f  U  +  n)  =  (c  -  Fz_)  -  (Fn)  (2.14) 

ij  5  c  •  Fz  ;  e^  =  F£  (2.15) 
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(2.3)  may  be  termed  as  the  'point  value-type'  quadrature  formula  as  it 
assumes  that  the  value  of  Ag  and  at  8j.  Xj  (or  perhaps  at  the  center  of  the  block 

ei  +  4r  «  A j  +  )  1S  constant  over  the  whole  block  a^..  We  may  easily  integrate 

Ynfn  over  the  block,  but  by  keeping  Zg  outside  the  integral  in  (2.5)  below  implies 
a  very  'smooth'  gravity  anomaly  signal  over  the  block  equal  to  the  mean  value. 

This  can  be  remedied  by  multiplying  by  a  'de-smoothing'  factor,  un  (Colombo, 

1981,  p.  33),  so  that  we  have  a  'mean  value-type'  quadrature  formula: 


where 


Un  N-12N-1  __ 

TTn^TT  fi0  jl0  4s1j  J°uVn>-  X) 


e  .+A6 


J  L(e.  x)da  =  I  Pnm(cos  0)  sine  de  / 


*j+AA  (cos 


x.  (sin 

J 


mx  dx  , 


(2.5) 


(2.6) 


The  evaluation  of  the  integral  of  associated  Legendre  functions  in  (2.6) 
will  be  discussed  in  Section  3.1.  The  de-smoothing  factor  Mn  may  be  related  to 
the  Pellinen/Meissl  smoothing,  or  averaging,  operator  en  (Rapp,  1977).  By 
numerical  tests  for  blocks  of  5°x  5°  and  larger,  Colombo  (1981,  p.  76,  (3.9))  found 
a  simple  nearly  optimal  relationship: 


pn  ~  4trnn  ’  nn 


6*  ,  0<n<N/3 
8  ,  N/3  <  n  <  N 


('  "n  ’ 

1  ,  n 


(2.7) 


The  use  of  various  de-smoothing  factors  in  (2.7)  may  cause  a  sharp  discon¬ 
tinuity  in  the  potential  coefficients  at  degree  N/3,  and  at  degree  N,  as  was 
noticed  by  Rapp  (1981,  p.  24).  Further,  pn  takes  into  account  only  the  block 
size,  as  en  is  given  by: 


Bn  *  cot 


3ni  (cos  »o) 
n(n+l) 


(2.8) 


where  it0  is  the  radius  of  a  circular  cap  having  the  same  area  aii  of  the  block. 

But  pn  does  not  take  into  account  the  error  estimate  of  Ag  employed  in  (2.5).  Hence, 
Colombo  (1981)  also  considered  the  rigorous  determination  of  'optimal  quadrature 
weights',  ,  based  on  the  minimization  of  the  sum  of  two  quadratic  terms,  one 
due  to  'sampTing  error',  and  the  other  due  to  'data  noise'.  This  is  discussed 
in  Section  2.2. 
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2.  HARMONIC  ANALYSIS  AND  SYNTHESIS 


We  will  henceforward  consider  l°x  1°  mean  anomalies,  referred  to  the  Geodetic 
Reference  System  1980  (Moritz,  1980b),  as  the  specific  gravity  related  quantity, 
as  in  (1.3).  We  will  also  adopt  a  more  compact  notation: 


c“ 

nm 


Yn>’ 


C  ,  as 

=0  / 

nm 

l 

S  ,  <*= 

=  1  ( 

nm 

Pnm(cos 

nm 

e) 

cos 

mx , 

a=0 

P_m(cos 

nm 

e) 

sin 

mx , 

a  =  l 

(2.1a) 


(2.1b) 


where  we  have  expressed  the  functions  in  (2.1b)  in  terms  of  co-latitude  e,  instead 
of  geocentric  latitude  £' .  We  will^  henceforward  assume  as  in  (2.1a)  that  the 
potential  coefficients  C^m,  i.e.  {Cnm,  Snm},  are  the  residual  potential  coefficients 
to  the  ellipsoidal  field  of  GRS  80,  and  we  will  dispense  with  the  notation 
of  superscript  asterisk  to  denote  such  residual  quantities. 

2.1  Quadrature  Formulas  with  De-smoothinq  Factors 

We  may  then  rewrite  (1.3)  for  harmonic  analysis  as  follows: 


To  implement  (2.2),  we  first  replace  the  integral  by  summation  over  the  anomaly 

blocks;  and  the  estimate  Cnm  is  now  affected  by  the  'sampling  error’  because  of 
the  finite  area,  of  the  actual  anomaly  blocks,  in  the  numerical  quadrature 
formula: 


,  N-l  2N-1 

ra  -  1  y  y 

nm  4^7  (n-l)  ^  ^ 


49  <V  \j>  Vnm<V  Xj>  41j 


(2.3) 


A.j j  =  ax  (cose^  -  cos(e..+Ae))  =  for  all  j  in  any  latitude  band  i, 

AX  =  Xj+.-Xj,  A6  =  e.+1-e.;  ag  =  ax  =  1°,  N  =  ~  =  180  (2.4) 


where  the  numerical  values  of  ag,  ax,  N  have  been  shown  for  l°x  1°  equi-angular 
blocks. 
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the  correct  order  of  kJim  required  is  the  vector  form  of  the  matrix  stored  column 
by  column.  This  would  then  allow  the  right  hand  side  vectors  to  be  read  in  suc¬ 
cessively  as  180  numbers  at  a  time  for  a  given  R(m)  matrix  (m=0  to  NN),  for  solving 
180  numbers  of  xnm  for  the  given  m,  for  degrees  n=m  to  NN. 

This  reordering  of  the  right  hand  side  vectors,  analogous  to  the  transpose 
of  a  matrix,  is  a  trivial  problem  if  the  whole  matrix  could  be  indexed  as  one 
array.  In  fact,  by  the  exploitation  of  symmetries  in  (3.3),  it  would  be  sufficient 
to  -reorder  the  right  hand  side  vectors,  obtained  from  I^m  in  the  manner  indicated 
above  for  i=0  to  89;  and  after  reordering  the  90  elements  of  each  right  hand  side 
vector,  extend  it  to  180  elements  using  (3.3).  However,  the  array  dimension 
required  for  reordering  the  right  hand  side  by  indexing  is  not  possible  for  the 
case  of  l°x  1°  anomalies,  as  is  apparent  from  Table  3.3. 


Table  3.3  Size  of  Array  Needed  for  Reordering  the  Right  Hand  Side  Vectors  by 
Indexing  for  Estimating  Coefficients  to  Degree  NN. 


Degree 

NN 

#  Elements  in 

Array  Dimension 

(Double  Precision) 

in  Mega  Bytes 

(NN+1) (NN+2)/2 

each  i 

Total 

i=0  to  89 

180 

16,471 

1,482,390 

11.86  MB 

250 

31,626 

2,846,340 

22.77  MB 

300 

45,451 

4,090,590 

32.72  MB 

The  problem  of  reordering  was  solved 
elements  of  were  correctly  ordered  in 
NN,  they  were  tagged  with  a  number  k. 


by  using  IBM  utility  SORT,  after  the 
a  latitude  band  i,  m=0  to  NN,  n=m  to 


k  =  j  *  100  +  (i+1) ;  j=l  to  NLL,  i=0  to  89;  NLL  =  (NN+l)(NN+2)/2  (3.7) 


and  written  on  tape.  After  the  tagged  numbers  were  sorted  in  ascending  sequence 
by  SORT,  they  were  in  the  order  needed  in  (2.30).  The  number  tag  was  stripped, 
and  the  right  hand  side  vector  _knm  written  out  on  another  tape.  This  step  of 
reordering  may  be  termed  as  step  3B.  The  CPU  time  by  the  SORT  routine  for 
sorting  2.8  million  elements  of  for  NN=250  on  an  Amdahl  470  V/8  computer 
was  2  min.  58  sec.,  while  the  initial  tagging  of  numbers  took  3  min.  2  sec.,  and 
the  writing  out  of  sorted  numbers  took  2  min.  33  sec. 

It  may  be  noted  that  steps  2,  3A  and  3B  do  not  depend  on  the  l°x  1°  anomaly 
error  estimates,  or  the  anomalies,  but  depend  on  the  degree  NN  up  to  which  co¬ 
efficients  are  to  be  estimated.  The  testing  of  the  effect  of  different  anomaly 
errors,  or  different  anomaly  data  sets  may  therefore  be  tried  for  low  values  of 
NN  like  12,  or  60,  or  180,  before  developing  coefficients  to  the  highest  degree. 
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3.33  Solution  of  the  Normals 

The  solution  of  the  normals  in  (2.30)  proceeds  by  reading  in  the  reformatted 
normals  R(m)  in  step  3A,  which  is  the  vector  form  of  the  upper  diagonal  of  symmetric  ^ 
180  x  180  matrix.  The  average  variance  o?  in  (2.32)  of  anomalies  in  each  latitude  * 

band  i,  i=0  to  179,  is  added  to  the  diagonal  elements  of  R(m),  yielding  the  matrix 
(R(m)+W)  in  (2.30).  This  is  done  for  each  R(m)  matrix,  one  at  a  time,  for  m=0  to 
NN . 

After  inversion  of  (R(m)+W),  90  numbers  of  the  reordered  right  hand  side  ' 

vectors  in  step  3B  are  read  in,  expanded  to  180  numbers  by  utilizing  the  symmetries 
in  (3.3),  and  the  optimal  quadrature  weights  Xnm  are  solved  for  180  numbers  at  a 
time,  x,L,  1=0  to  179.  The  solution  proceeds  for  each  m=0  to  NN,  n=m  to  NN,  re¬ 
quiring  (NN+l)  matrix  inversions,  and  (NN+l)(NN+2)/2  solutions.  The  CPU  time 
for  this  step  increases  sharply  for  large  NN.  An  Amdahl  470  V/8  computer  took 
about  34  minutes  CPU  time  for  this  step  at  NN=250.  » 

The  total  error  variance  estimates  in  (2.21)  are  computed  by: 


2_  _  a2_  _  1  cn 

°  cnm  snm  x2(n-l)2  2n+l 


AA' 


1-COSJP'AA  .f 
_^2  '  ‘ 


m^ 


if  m=0 
ntfO 


N-l 
2N  ■  l 


L  k1  X1 
i=0  nm  nm 


(3.8) 


which  follows  from  the  similarity  of  terms  in  (2.20)  and  (2.21);  and  the  expres¬ 
sions  (2.31)  and  (2.26)  using  x^  and  k^.  The  variances  were  earlier  computed 
only  for  a  pair  of  coefficients,  unless  m=0,  as  discussed  below  (2.32).  A  modi¬ 
fication  was  made  to  compute  the  standard  error  of  each  coefficient;  and  the 
factor  of  h  for  ntfO,  for  making  the  variance  for  each  coefficient  in  a  pair  equal, 
i.e.  a2.nm  =  has  been  incorporated  in  (3.8) 

The  error  variance  [Epltrnm  =  t^n]-5-nm,  for  each  coefficient  due  to  propagated 
noise  is  computed  similarly  to  (3.8)  by: 


^n-^Cnm  ^n^snm 


2n+l 


AA2  if  m=0 

hsmm  if  m 

m  2 


N-l 


•  2N  •  l  x! 


.2 


i=0 


nm 


(3.9) 


And,  the  error  variances  due  to  sampling  finite  size  of  l°x  1°  anomaly  blocks  is 
computed  by: 


^Es-^rnm 


^Es^nm 


=  a2 


rnm 


^n^nrn 


-  a. 


Snm 


^n-^nm 


(3.10) 


The  program  sums  up  the  variances  per  degree  due  to  sampling,  propagated 
noise,  and  total.  Modification  was  made  to  also  compute  other  statistics  for 
corresponding  standard  errors  per  degree,  percentage  sampling  and  noise  errors 
of  the  total  error,  and  also  percentages  of  the  errors  per  potential  coefficient 
variation,  on,  per  degree: 


I 


l 


i 
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(3.11) 


on  =Vc^>(  y(n-l)) 


This  step  of  solving  the  normals  for  computing  the  optimal  quadrature  weights 
and  the  error  estimates  for  the  coefficients,  based  on  average  anomaly  variances 
in  different  latitude  bands,  may  be  termed  as  step  3C. 

3.4  Computation  of  Potential  Coefficients 

The  computations  of  potential  coefficients  in  (2.31),  from  the  optimal  quadra 
ture  weights  iLnm  obtained  in  step  3C,  is  done  by  the  harmonic  analysis  of  global 
anomaly  data  by  subroutine  ANALYS.  The  computations  are  shortened  by  recognizing 
that  the  symmetries  in  (3.3)  also  apply  to  ^,m.  Hence,  (2.31)  is  transformed  to: 


The  implementation  of  (3.12)  is  further  speeded  up  by  computing  ai,  b^,  a{Jj"  "  , 
bjj}-^  as  the  Fourier  transform  along  latitude  bands  i  and  N-l-i  of  the  values 
of  Ag  (e,  a)  ,  i .e. 


fcosj 
sin) 


mj  aa  Ag  (ei ,  Aj) 


2N-1 

l 

j=0 


mj  aa  Ag  (e^j.-j.  \j) 


(3.13) 


A  modification  was  made  to  ANALYS  to  compute  several  sets  of  coefficients 
in  the  same  run,  either  combining  several  sets  of  optimal  quadrature  weights 
associated  with  different  anomaly  error  estimates  with  one  global  anomaly  data 
set,  or  by  utilizing  different  anomaly  data  sets  with  one  set  associated 
with  a  certain  set  of  anomaly  errors.  The  CPU  time  for  the  harmonic  analysis 
of  64,800  l°x  1°  anomalies  to  estimate  two  potential  coefficient  sets  up  to  de¬ 
gree  NN=250  was  only  1  min.  26  sec.  on  an  Amdahl  470  V/8  computer.  This  step 
of  computing  the  potential  coefficients  may  be  termed  as  step  4. 


There  were  two  additional  steps  of  the  comparison  of  potential  coefficient 
sets  by  degree,  and  of  computing  error  estimates  by  degree  and  cumulatively  for 
the  undulations  and  anomalies  implied  by  a  coefficient  set.  These  may  be  termed 
as  steps  5  and  6,  and  are  discussed  in  Sections  5  and  6  respectively. 

The  various  steps  in  the  algorithm  for  optimal  estimation  of  coefficients 
are  summarized  in  Table  3.4. 
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4.  POTENTIAL  COEFFICIENT  AND  GRAVITY  ANOMALY  DATA  SETS 


The  parameters  affecting  the  optimal  estimation  of  coefficients  are  given 
in  the  last  column  of  Table  3.4.  The  highest  degree  Nmax,  to  which  the  summations 
are  carried  out  for  computing  the  Fourier  transforms  of  the  covariances  of  l°x  1° 
anomalies  in  (3.2),  was  kept  as  300  based  on  experiments  carried  out  by  Colombo 
(1981,  Sec.  4.2).  The  global  anomaly  degree  variance  model  used  by  Colombo  (ibid.. 
Sec.  3.1)  was  also  retained.  It  consisted  of  empirically  obtained  anomaly  degree 
variances  up  to  degree  100  from  a  180  x  180  potential  coefficients  field  developed 
by  Rapp  in  1978,  and  by  a  two  component  variance  model  described  by  Rapp  (1979, 
p.  15,  Table  5,  Case  1)  for  degree  n >  100: 


(4.1) 


o  =  3.4050,  a ^  =  140.03,  s  =  0.998006,  s  =  0.914232,  A  =  1,  B  =  2 

The  other  parameters,  and  data  sets,  used  in  the  optimal  estimation  of  coefficients 
are  described  below. 

4.1  Anomaly  Error  Estimates 


Several  different  sets  of  anomaly  error  estimates  were  used.  One  set  computed 
average  variance,  a?  in  (2.32),  in  each  latitude  band  from  error  estimates  of  a 
global  l°x  1°  anomaly  data  set  utilizing  terrestrial  anomalies  merged  with  anomalies 
determined  from  Seasat  data,  and  balance  anomalies  computed  from  potential  coef¬ 
ficients  being  assigned  standard  error  of  30  mgals.  This  anomaly  data  set  was  a 
slightly  updated  version,  with  some  additional  terrestrial  data,  of  the  set  described 
in  Section  2.3.  The  computed  values  of  average  anomaly  standard  errors,  on- ,  per 
latitude  band  are  given  in  Table  4.1,  and  this  will  be  termed  as  anomaly  error 
estimate  A.  The  standard  errors  range  from  30  to  about  4  mgals. 


Table  4.1  Average  Anomaly  Errors  in  1°  Latitude  Bands  for  Error  Estimate  A. 


Latitude 

Bands 

i 

Average  Anomaly  Standard  Error, 

in 

mgals 

0 

to 

14 

30 

27 

27 

27 

27 

27 

28 

27 

27 

26 

25 

25 

24 

23 

22 

15 

to 

29 

22 

21 

20 

19 

19 

16 

16 

15 

16 

16 

15 

15 

15 

15 

14 

30 

to 

44 

12 

12 

11 

11 

11 

11 

11 

11 

10 

11 

10 

9.9 

10 

10 

11 

45 

to 

59 

11 

11 

11 

11 

12 

12 

11 

11 

12 

12 

10 

10 

10 

11 

10 

60 

to 

74 

11 

10 

9.8 

10 

9.9 

10 

10 

9.2 

9.4 

9.7 

9.7 

8.9 

9.3 

8.3 

8.1 

75 

to 

89 

6.9 

7.1 

8.3 

8.5 

8.9 

9.1 

8.8 

8.8 

8.8 

9.1 

8.8 

7.5 

7.6 

7.9 

7.4 

90 

to 

104 

8.0 

8.2 

8.4 

8.4 

8.9 

9.1 

8.7 

8.3 

8.7 

7.7 

7.3 

7.3 

7.7 

7.4 

7.6 

105 

to 

119 

8.1 

7.4 

7.4 

7.8 

7.2 

7.6 

7.2 

6.8 

6.9 

6.5 

6.4 

6.0 

5.6 

5.7 

5.8 

120 

to 

134 

5.6 

5.3 

5.0 

4.2 

4.6 

3.8 

4. 1 

3.9 

4.6 

4.9 

4.1 

4.1 

4.4 

4.8 

5.2 

135 

to 

149 

4.4 

4.0 

3.6 

3.9 

4.4 

5.1 

5.0 

4.4 

4.5 

4.6 

3.8 

3.6 

3.3 

3.7 

5.5 

150 

to 

164 

3.7 

3.9 

4.5 

5.2 

5.9 

6.7 

8.0 

11 

12 

13 

27 

28 

27 

28 

28 

165 

to 

179 

26 

26 

25 

24 

24 

24 

24 

24 

23 

25 

26 

26 

23 

14 

14 

A  second  set,  termed  anomaly  error  estimate  B,  assigned  uniformly  for  each 
1°  latitude  band  average  anomaly  standard  error,  aj,  of  5  mgals.  Error  estimate 
A  is  a  realistic  estimate  of  the  current  averaged  anomaly  variances  in  1°  latitude 
bands  for  l°x  1°  anomalies.  Error  estimate  B  assumes  no  latitude  variation  in 
averaged  variances,  and  is  therefore  an  idealized  version  of  anomaly  error  distri¬ 
bution.  A  value  of  5  mgals  is  an  optimistic  estimate  giving  a  lower  error  bound 
till  the  Geopotential  Research  Mission  is  flown  at  the  end  of  this  decade. 

A  uniform  average  anomaly  standard  error  of  20  mgals  was  also  used  in  some 
tests,  as  a  20  mgal  error  was  used  for  each  anomaly  in  generating  the  current 
Rapp  (1981)  coefficients,  as  discussed  at  the  end  of  Section  2.3  of  this  study. 
Another  set  of  error  estimates  was  developed  on  the  same  lines  as  error  estimate 
A,  but  only  using  terrestrial  data  resulting  in  a  range  of  average  standard  error, 
ai,  from  30  to  9  mgals.  Some  other  anomaly  error  estimates  would  be  described 
in  the  tests  in  Sections  5.1  and  5.2. 

4.2  Data  Sets  A 

The  development  of  the  current  set  of  merged  potential  coefficients  to  degree 
and  order  180  (Rapp,  1981)  was  described  in  Section  2.3.  This  will  be  termed  as 
coefficient  set  A,  and  will  be  compared  with  other  coefficient  sets  developed 
by  optimal  estimation. 

Coefficient  set  A  was  used  to  compute  a  global  set  of  64,800  l°x  1°  mean 
gravity  anomalies  by  (2.11)  using  subroutine  SSYNTH,  with  y=979, 800  mgals.  This 
global  set  will  be  termed  as  anomaly  set  A.  As  the  maximum  degree,  Nmax,  in 
the  summation  in  (2.11)  was  180,  the  anomaly  set  A  can  only  be  used  for  estimating 
coefficients  to  the  highest  degree  NN  of  180.  But  the  anomaly  set  A  would  be 
a  good  data  set  to  test  the  effect  of  anomaly  error  estimates  on  the  optimal 
estimation  of  coefficients,  by  testing  the  latter  against  coefficient  set  A. 

These  tests  will  be  described  in  Section  5.1. 

4.3  Data  Sets  B 

Anomaly  set  B  was  the  global  set  of  64,800  l°x  1°  'adjusted  anomalies'  described 
in  Section  2.3  obtained  by  implementing  (2.38).  Anomaly  set  B  therefore  retains 
the  high  degree  (n  >  180)  information  present  in  the  terrestrial  and  the  Seasat 
altimeter  derived  anomalies,  which  were  used  to  develop  the  'adjusted  anomalies'. 

As  mentioned  in  Section  2.3,  the  coefficient  set  developed  by  the  harmonic 
analysis  of  anomaly  set  B,  by  implementing  (2.5)  with  de-smoothing  factors  (2.7), 
using  subroutine  HARMIN  will  be  termed  as  coefficient  set  B.  Coefficient  set  B 
would  be  considered  only  up  to  degree  and  order  180. 

Coefficient  sets  A  and  B  are  the  same  from  degree  49  to  180.  'Adjusted 
SET  1 ’  coefficients,  complete  to  degree  and  order  36,  with  some  additional  coef¬ 
ficients  to  degree  48,  replaced  the  corresponding  coefficients  in  coefficient 
set  B  to  form  the  coefficient  set  A.  The  difference  between  these  two  sets  is 
shown  in  Table  5.1. 

The  coefficients  developed  by  optimal  estimation  from  anomaly  set  B  are  there¬ 
fore  comparable  to  coefficient  set  B  through  degrees  2  to  180,  or  to  coefficient 
set  A  through  degrees  49  to  180.  The  comparison  with  coefficient  set  A  through 
degrees  2  to  48  would  also  show  the  additional,  though  small,  effect  due  to  the 
simplified  adjustment  model  in  (2.35). 


5.  COMPARISON  OF  SETS  OF  POTENTIAL  COEFFICIENTS 


If  the  potential  coefficients  in  two  different  estimates  to  degree  NN  are 
denoted  by  (Cpm,  Snm)  and  (C^m,  S^m),  then  the  two  sets  of  estimates  may  be  com¬ 
pared  in  magnitude  by  percentage  coefficient  difference  (%An)  by  degree  n,  undu¬ 
lation  difference  ( /un)  by  degree,  and  anomaly  difference  (/Acn)  by  degree. 


Aon 

%A  -  — — 
n  °n 


100 


(5.1) 


n 


o‘  *  y  (c2  +  s2  ) 

L  -  nm  nm' 


i  ■ 


m=0 

n 


Ao2  =  l  [(C  -  C')2  +  (S 

n  _  «  nm  nm 
m-0 


nm 


S*  )2] 
nm'  J 


(5.2) 

(5.3) 


/a i  =  R  •  A a  meters  (5.4) 

n  n  v  ' 

-  y(n-l)  •  Aan  mgals  (5.5)  •;!; 


where,  as  before,  y  =  979,800  mgals  is  the  average  value  of  gravity  over  the  whole  |> 

globe  approximated  by  a  sphere  of  radius  R=6.371xl06  meters.  For  example,  co- 
efficient  sets  A  and  B  differ  up  to  degree  36;  differ  in  some  additional  coefficients 
up  to  degree  48,  and  are  the  same  in  higher  degrees  to  180.  This  comparison  is  <<• 

shown  in  Table  5.1.  V-V^ 

Table  5.1  Comparison  of  Potential  Coefficient  Sets  A  and  B  (Sec.  4.2  and  4.3), 


Degree  n 

%  Coefficient 
Difference 

Undulation 

Difference 

in  meters 
n 

Anomaly 
Difference 
/ac^  in  mgals 

1 

2 

1.23 

.22 

.03 

3 

.34 

.06 

• 

4 

.39 

.04 

.02 

5 

.64 

.05 

.03 

6 

.51 

.03 

.02 

9 

.54 

.01 

.02 

12 

.56 

.01 

.01 

\V.\« 

24 

.51 

.00 

.01 

36 

.71 

.00 

.01 

37 

.03 

.00 

.00 

V  V 

42 

.17 

.00 

.00 

48 

.06 

.00 

.00 

.-v-v 

49-180 

.00 

.00 

.00 

•  .i 

5.1  Initial  Tests  with  Data  Set  A 


We  first  show  in  the  first  set  of  columns  in  Table  5.2,  a  comparison  between 
coefficient  set  A  and  another  coefficient  set  computed  to  degree  180  by  harmonic 
analysis  of  the  anomaly  set  A  by  HARMIN  in  (2.5),  using  de-smoothing  factors  in 
(2.7).  The  differences  would  arise  primarily  due  to  sampling  error  because  of 
using  l°x  1°  anomaly  blocks  of  a  finite  size,  instead  of  using  an  infinite  number 
of  point  anomalies  implied  by  the  integral  in  (2.2).  Additional  differences 
would  also  arise  due  to  any  inadequacy  in  the  de-smoothing  factors,  particularly 
a  discontinuity  in  (2.7)  near  degree  N/3=60,  which  is  seen  in  Table  5.2.  We 
also  see  an  increase  in  the  sampling  error  with  higher  degrees.  Anomaly  error 
estimates  are,  of  course,  not  considered  in  (2.5). 

In  the  next  two  sets  of  columns  in  Table  5.2,  coefficient  set  A  is  compared 
with  coefficient  sets  developed  from  anomaly  set  A  by  optimal  estimation  procedure, 
using  error  estimate  B  (a  =5  mgals),  and  error  estimate  A  (a.  ranging  from  30 
to  4  mgals  in  Table  4.1)  Respectively.  We  would  now  see,  besides  the  effect  of 
sampling  error,  larger  differences  because  of  propagated  noise  of  the  anomalies. 

We  would  also  expect  that  the  larger  the  uncertainty  in  anomaly  error  estimate, 
the  larger  would  be  the  difference  from  coefficient  set  A.  We  do,  in  fact,  see 
increased  differences  in  the  second  set  of  columns  over  the  first  set  of  columns, 
and  further  increased  differences  in  the  third  set  of  columns  over  the  second 
set.  However,  we  do  not  find  any  discontinuity  in  the  optimal  estimation  procedure 
in  the  second  and  third  sets  of  columns  near  degree  60,  as  is  seen  in  the  first 
set  of  columns  with  de-smoothing  factors. 

Table  5.2  Variation  of  Coefficients  from  De-smoothing  Factors  with  Coefficients 
from  Optimal  Estimation. 

Comparison  of  Coefficient  Set  A  with  3  Coefficient  Sets  Obtained  from 
Anomaly  Set  A. 


De-smoothing  Factors 
(2.7)  with  (2.5) 

Optimal  Estimation  Procedure 
(Table  3.4) 

Degree 

Anomaly  Errors 

Anomaly  Error 

Anomaly  Error 

n 

Not  Considered 

Directly 

Estimate  B** 

Estimate 

A** 

%An 

/AZn  m 

mgal 

%An 

Af  n  m  Acn  mgal 

%An 

m 

/Acn  mgal 

2 

.61 

.11 

.02 

.64 

.11  .02 

.85 

.15 

.02 

3 

.44 

.08 

.03 

.48 

.09  .03 

.40 

.08 

.02 

6 

.56 

.03 

.02 

.58 

.03  .03 

.63 

.04 

.03 

12 

.91 

.01 

.02 

.87 

.01  .02 

1.4 

.02 

.03 

36 

1.2 

.00 

.02 

1.7 

.00  .02 

5.0 

.01 

.07 

60 

2.1 

.00 

.04 

2.3 

.00  .04 

6.9 

.01 

.13 

61 

4.4 

.01 

.07 

2.9 

.01  .05 

8.6 

.02 

.14 

90 

8.9 

.01 

.14 

6.0 

.01  .10 

15.2 

.02 

.25 

120 

14.2 

.01 

.23 

10.7 

.01  .17 

21.4 

.02 

.34 

150 

22.3 

.01 

.30 

22.2 

.01  .30 

33.9 

.02 

.45 

180 

27.4 

.01 

.35 

36.2 

.02  .46 

46.8 

.02 

.60 

^Anomaly  errors  were  not  considered  in  computing  coefficients  from  anomalies  using 
de-smoothing  factors. 


♦♦Error  Estimate  B:  o1  in  (2.32)  =  5  mgals. 

Error  Estimate  A:  o,  ranges  from  30  to  4  mgals,  see  Table  4.1. 
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Table  5.3  Variation  of  Optimally  Estimated  Coefficients  with  Variation  in  Anomaly  Error  Estimates. 
Comparison  of  Coefficient  Set  A  with  6  Coefficient  Sets  Estimated  from  Anomaly  Set  A. 


- 1 1<? 
ic  l< 
01  S. 

VO  E  - 

0)0  lo* 

l/l  C\J  [< 


ID 

31  2 

■5  8 


C  r- 
U  <V 


Ol  I1? 


AS 

*  01 

co  E 


"  c 

< 

o  « 


C  r- 

O  <o 
_  <a  cn 
'Z  S.  E 


C\J  g1  - 

E  I  c 

ai  _ 

w>  U3 

<0  .1  > 


-  1“ 

ra  L< 

cn  S. 

e  _ 

fH 

—v  I  C 
qj  O 

vn  O  l< 
ns  • 


coronmunvor'-r-.o 

OOOOO 


OOVOOOVOtOLO^S-vO 

CVJ^HOOOOOOO 


CVJ  VO 
LO  VO 


00000 


l£> 

cn 

cn 

fH 

in 

10 

CM 

CM 

LO 

<tcomsinmc\jwn 

■— 100000000 


OOvfrteOOONO^H 

N<3-lflO)CDO)'-iHS 


CMCVjCVJ^J-COCNjCOCVJfO 

OOOOOOOOO 


inoO'JiaVNNHN 

OOOOOOOOO 


mo^Mnixo<'J'5’ 

ajvtvtooiO'Ofl'N^ 


CMnM^-cnMMevjcvj 

OOOOOOOOO 


•-HCD^-f^.cOCMCM'-v^H 

OOOOOOOOO 


*a-oovoroooo>r^oor^. 


CMCOCVJ^-CNJfOCVJCJH 

OOOOOOOOO 


rtSl^tNMNHHH 

OOOOOOOOO 


iH  00  vf>  CVJ  VO  CvJ  00  co 

(DvfvfWlflNlOlO® 


PJfON<rNCVJN(*)H 

OOOOO 


«-t  00  vo  co  1  1— 1 

1-1  O  O  O  O  )  o 


f — 1  lt>  in  . — — 1  cvj  n 

vD<J-^J-00Lf)r^vOOC0 


m  in  in  vo  O) 


OHM 


To  examine,  in  greater  detail,  the  effect  of  anomaly  error  estimates  on 
coefficients  by  optimal  estimation,  several  different  sets  of  coefficients  were 
computed  with  anomaly  error  estimates  as  shown  in  Table  5.3,  which  also  shows 
the  comparison  with  coefficient  set  A.  There  is  some  instability  in  the  inverse 
of  matrix  (R(m)+W)  in  (2.30),  when  the  elements  of  diagonal  matrix  W  in  (2.32) 
correspond  to  aj  of  .001  mgals.  But  with  °-j21  mgal ,  Table  5.3  shows  conclusively 
that  optimally  estimated  coefficients  depend  on  the  magnitude  of  anomaly  error 
estimates. 

The  tests  in  Table  5.3  were  run  with  NN=12,  after  ensuring  that  the  estimated 
coefficients  are  not  correlated,  and  hence  the  results  in  Table  5.3  with  NN=12 
are  also  valid  for  higher  degrees.  This  conclusion  was  arrived  at  by  repeating 
the  optimal  estimation  of  coefficients  for  anomaly  set  A,  with  both  error  estimates 
A  and  B,  for  NN=180  in  Table  5.2,  also  for  NN=60  and  NN= 12 .  The  common  coef¬ 
ficients  were  the  same  whether  the  coefficients  were  estimated  to  NN= 12 ,  60,  or 
180.  Similar  results  were  obtained  later  with  anomaly  set  B. 

5.2  Initial  Tests  with  Data  Set  B 

Tests  were  first  made  to  examine  the  variation  in  optimally  estimated  coef¬ 
ficients  from  anomaly  set  B  by  using  different  error  estimates  as  in  Table 
5.3,  and  comparing  the  optimally  estimated  coefficients  against  the  coefficient 
set  B.  The  results  were  similar  to  those  presented  in  Table  5.3  showing  larger 
variation  in  the  optimally  estimated  coefficients  as  larger  anomaly  error  esti¬ 
mates  were  used,  reflecting  the  effect  of  propagated  anomaly  errors. 

It  was  also  noticed  that  generally  there  was  a  better  agreement  between 
coefficient  set  B  and  optimally  estimated  coefficients  from  anomaly  set  B  for 
any  given  anomaly  error  estimate,  as  contrasted  with  the  case  of  agreement  between 
coefficient  set  A  and  optimally  estimated  coefficients  from  anomaly  set  A  pre¬ 
sented  in  Section  5.1.  For  example,  this  may  be  seen  in  Table  5.4  for  the  anomaly 
error  estimates  for  a-j=  5  mgals,  and  1  mgal ,  when  contrasted  with  similar 
results  in  Table  5.3.  This  may  be  due  to  the  anomalies  in  set  B  having  been 
obtained  in  a  combined  adjustment  with  ' SET 1  *  coefficients  (Section  2.3),  yielding 
the  anomaly  set  B  and  the  'adjusted  SET1 '  coefficients.  The  adjustment  of  anomalies 
in  set  B  not  only  adjusts  the  anomaly  spectrum  to  degree  36  (including  some  co¬ 
efficients  to  48),  but  perhaps  also  modifies  to  some  extent  the  higher  degree 
spectrum.  This  will  be  seen  in  the  results  in  Section  5.4. 

However,  we  also  find  curiously  a  larger  disagreement  in  the  comparisons 
in  Section  5.2  at  the  very  low  degrees,  when  contrasted  with  the  comparisons 
in  Section  5.1.  This  may  be  seen  again  in  Table  5.4  for  the  anomaly  error  esti¬ 
mates  for  o-j  =  5  mgals,  and  oi  =  1  mgal,  when  contrasted  with  similar  results 
in  Table  5.3  for  degrees  2  and  3.  This  may  be  due  to  the  simplified  adjustment 
model  in  (2.35)  and  (2.36)  being  generally  adequate  at  degrees  greater  than  3. 

Or,  this  may  be  an  artifact  of  forcing  the  weight  matrix  of  anomalies  in  (2.37) 
as  cos  <j )'/o2,  while  the  weight  matrix  in  optimal  estimation  is  implicitly 
1/a?  in  (2.30),  read  with  (2.22). 

Several  sets  of  coefficients  were  then  optimally  estimated  with  oi  being 
modified  from  1  mgal  to  1/cos  <j>j  mgal,  and  l//cos  mgal  to  match  the  anomaly 
weights  in  (2.37)  as  cos  <t>'/a2.  However,  this  made  the  agreement  with  coefficient 
set  B  worse  for  all  degrees  as  compared  to  the  case  of  =  1  mgal.  Other  sets 
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of  coefficients  were  then  optimally  estimated  with  o-j  =  1  •  cos<j>.-  mgal ,  and 
1  •  fcos^y- mgal .  The  last  case  resulted  in  close  agreement  with  coefficient 
set  B  at  degrees  greater  than  3,  but  the  disagreement  persisted  at  degrees  2 
and  3,  and  was  now  of  a  larger  amount.  All  these  tests  were  run  with  coef¬ 
ficients  estimated  to  degree  NN=12. 

Table  5.4  shows  the  comparison  of  coefficient  set  B  with  four  optimally 
estimated  coefficient  sets  by  harmonic  analysis  of  anomaly  set  B.  The  anomaly 
error  estimates  in  the  four  cases  were  oj  =  5,  1,  1-  cos  <j>)  and  1  •  /cos  mgal. 
The  first  two  cases  are  for  comparison  with  similar  cases  in  Table  5.3  for  the 
data  set  A, and  the  last  two  cases  in  Table  5.4  were  an  attempt  to  reduce  the 
disagreements  at  degrees  2  and  3.  The  tests  were  repeated  with  data  sets  A  in 
Section  5.1  with  =  1/cos  <f>.| ,  1/  /cos  $  j ,  1,  1  •  cos  <f>j ,  and  1  •  /cos  <jT[  but 
did  not  show  any  larger  differences  as  noticed  in  the  case  of  data  sets  B. 

As  already  mentioned,  the  optimally  estimated  coefficients  from  the  adjusted 
anomalies  in  set  B,  adjusted  through  the  simplified  model  in  (2.35),  show  generally 
a  better  agreement  with  coefficient  set  B  as  contrasted  with  the  agreements  of 
data  sets  A  in  Section  5.1,  except  for  degree  2  and  3.  Coefficient  sets  will  be 
optimally  estimated  from  anomaly  set  B  to  higher  degrees  up  to  NN=25Q  in  Section 
5.4  for  the  cases  of  anomaly  error  estimates  A  and  B,  without  forcing  these  error 
estimates  to  ai  /cos  <j>{  or  a\/  /cos  <j>^ . 

5.3  Test  with  Altered  Anomaly  Data  Set 

The  assembling  of  64,800  l°x  1°  global  anomalies  was  described  in  Section 
2.3,  where  8049  anomalies  were  assigned  values  implied  by  1 SET1 ’  coefficients 
(complete  to  degree  36),  and  a  standard  error  of  30  mgals.  During  the  combined 
adjustment  of  anomalies  and  SET1  coefficients,  anomalies  with  higher  standard 
errors  received  larger  corrections.  The  location  of  3827  anomalies,  which  had 
corrections  larger  than  7  mgals,  were  shown  in  Rapp  (1981,  p.  25,  Fig.  6)  and 
is  now  reproduced  here  as  Figure  5.1. 

Four  10°x  10°  blocks  are  marked  in  Figure  5.1,  two  each  in  the  northern 
hemisphere  in  Central  Siberia  and  Central  Africa,  and  two  each  in  the  southern 
hemisphere  in  Southwest  Africa  and  Central  Andes  in  South  America.  The  lati¬ 
tudinal  and  longitudinal  limits  of  these  blocks  are  listed  in  Table  5.6.  Three 
tests  were  made  to  examine  the  effect  on  optimally  estimated  coefficients  if 
some  l°x  1°  anomal ies  in  the  global  anomaly  data  are  set  to  zero  instead  of  the 
value  implied  by  a-priori  coefficients.  From  the  global  anomaly  set  A,  three 
other  anomaly  sets  were  obtained:  anomaly  set  1  with  200  l°x  1°  anomalies  set 
to  zero  in  the  two  10°x  10°  blocks  in  the  northern  hemisphere;  anomaly  set  2  with 
200  1°X  1°  anomalies  set  to  zero  in  the  two  10°x  10°  blocks  in  the  southern 
hemisphere;  and  anomaly  set  3  with  400  l°x  1°  anomalies  set  to  zero  in  the  four 
10°x  10°  blocks  in  both  hemispheres.  Four  set:  of  coefficients  were  optimally 
estimated,  all  to  degree  NN=60,  from  the  global  anomaly  sets  A,  1,  2,  and  3. 

The  realistic  anomaly  error  estimate  A,  o-j  =  30  to  4  mgals  in  Table  4.1,  was 
used  in  all  four  cases.  The  comparison  of  coefficient  set  A  with  these  4  sets 
of  coefficients  is  shown  in  Table  5.5. 


Coefficient  Percentage  Error  Coefficient  Set  A 

As  A  Function  Of  Degree 
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6.2  Error  Estimates  by  Degree 

The  expression  for  error  variances  of  a  coefficient,  Cpm,  due  to  sampling 

error,  [E<r]_  ,  propagated  noise,  [En]_  m  ,  and  the  total  variance,  o2  ,  were 

b  Cnma  3  "  HTnma  nma 

given  in  (3.10),  (3.9)  and  (3.8)  respectively.  We  may  now  define  the  error  vari¬ 
ances  per  degree  due  to  sampling,  a2sn,  due  to  propagated  noise, a|nn,  and  the 
total  error  variance  o2n  per  degree,  as: 


f  °esn> 

l  i 

1  1 

i  fE-l  \ 

L  S-Tnma  | 

n  1 

j  / 

\°enn  | 

>=  I  I  < 

|  ^n^Cnma  ^ 

,  m^0  a=0  , 

°en ' 

)  \ 

where  we  have  used  [Et]^.  _  =  o2  in  view  of  (2.16),  and  [•]  indicates  one  element 

L  1  -xnma  nma  : 

of  the  matrix.  The  errors  (square  root  of  variance)  per  degree,  a  a  .  orrl, 

n  3  esn  cnn  en 

may  be  expressed  as  percentages  of  potential  coefficient  variation,  on,  per  degree 

in  (3.11).  If  we  denote  these  percentage  errors  due  to  sampling,  noise  and  total 

as  %S  ,  %N  ,  %Tn  per  degree,  then: 


%Sn| 

%N 

1  iaesn  j 

>  =  < cr  „  > 

t — ‘ 

o 

o 

(6.2) 

"  1 

\  enn  / 

%T 

n  / 

1  (aen  ) 

with  an  =  </c^/(Y(n-l)) 


(6.3) 


The  percentage  errors,  %Sn,  %Nn,  %Tn,  have  been  tabulated  in  Table  j.3  for 
coefficient  sets  B1  and  B2.  We  note  that  %Sn  is  small  at  low  degrees,  but  in¬ 
creases  sharply  with  degree.  The  sampling  error  predominates  after  around  degree 
120,  is  more  than  50%  of  the  coefficient  value  at  degree  180,  and  more  than  80% 
of  the  coefficient  value  at  degree  250.  Because  of  the  predominant  effect  of 
the  sampling  error  at  high  degrees,  we  find  that  the  total  percentage  error, 

%Tn,  is  only  slightly  different  for  sets  B1  and  B2  at  degree  250,  though  it  is 
much  smaller  at  lower  degrees  for  set  B2.  However,  %Tn  does  not  reach  100%  in 
any  case  even  at  degree  250,  because  the  optimal  estimator  cannot  have  a  larger 
error  than  a  null  estimator  (which  will  correspond  to  100%  error),  (Colombo, 
1981,  p.  73). 
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Figure  6.2 
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Degree 


Table  6.2  Undulation  and  Anomaly  Magnitudes  per  Degree  and  Cumulatively. 

Comparison  of  New  and  Current  Sets  of  Potential  Coefficients. 
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Coefficient  Set  B1  estimated  optimally  from  anomaly  set  B  to  degree  250  with  anomaly  error  estimate  A  : 
ranging  from  30  to  4  mgals  (see  Table  4.1). 

Coefficient  set  A  is  the  current  set  of  coefficients  (Rapp,  1981)  to  degree  180.  (Also  see  Sec.  2.3). 
Anomaly  errors  were  not  considered  in  computing  coefficients  from  anomalies  using  de-smoothing  factors. 
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Table  6.1  Variation  of  Optimally  Estimated  Coefficients  to  Degree  250  with 
Variation  in  Anomaly  Error  Estimate. 

Comparison  of  Coefficient  Set  B1  with  Coefficient  Set  B2. 


Degree 

n 

%  Coefficient 
Difference 
%An 

Undulation 

Difference 

/A£n  in  meters 

Anomaly 
Difference 
/Acn  in  mgals 

2 

3.4 

.61 

.09 

3 

9.2 

1.70 

.52 

6 

4.0 

.23 

.17 

12 

10.5 

.11 

.18 

24 

6.6 

.03 

.09 

36 

4.4 

.01 

.06 

60 

5.5 

.01 

.09 

90 

11.7 

.01 

.18 

120 

14.6 

.01 

.22 

150 

20.2 

.01 

.23 

180 

26.8 

.01 

.24 

210 

30.6 

.01 

.22 

240 

37.0 

.01 

.20 

250 

40.6 

.00 

.19 

Both  coefficient  sets  B1  and  B2  were  estimated  from  anomaly  set  B  to  degree  and 
order  250  using  anomaly  error  estimates  A  and  B  respectively. 

Error  Estimate  A  :  ai  ranges  from  30  to  4  mgals  (see  Table  4.1). 

Error  Estimate  B  :  a.  =  5  mgals. 
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6.  OPTIMAL  ESTIMATION  OF  COEFFICIENTS  TO  DEGREE  AND  ORDER  250 


As  the  anomaly  set  A  was  computed  by  (2.11)  from  coefficient  set  A  to  degree 
Nmax=180,  t  could  not  be  used  to  compute  coefficients  beyond  degree  180.  Accordingly 
anomaly  set  B  was  used  to  optimally  estimate  coefficients  to  degree  NN=2 50,  using 
both  anomaly  error  estimates  A  and  B.  We  will  term  these  coefficient  sets  B1  and 
B2  respectively.  As  the  anomaly  error  estimate  A  represents  our  current  knowledge 
of  l°x  1°  anomalies  realistically,  coefficient  set  B1  would  represent  the  current 
estimates  of  high  degree  gravity  field.  However,  comparisons  with  coefficient  set 
B2  would  give  us  the  upper  bound  of  improvement,  with  respect  to  anomaly  error  esti¬ 
mates,  that  we  may  expect  until  the  end  of  this  decade. 

6.1  Magnitude  Information  by  Degree 

We  first  list  in  Table  6.1  the  percentage  coefficient  difference,  undulation 
difference,  and  anomaly  difference  by  degree,  £An,  /A2n ,  Acn,  between  coefficient 
sets  B1  and  B2.  We  notice  large  differences  due  to  different  anomaly  error  esti¬ 
mates  A  and  B,  as  we  would  expect  from  the  information  to  degree  180  in  Table 
5.7.  The  disagreement  at  degree  3  is  particularly  noticeable  because  of  strong 
equatorial  assymmetry  of  error  estimate  A  in  Table  4.1. 

We  next  list  in  Table  6.2,  the  undulation  magnitude  and  anomaly  magnitude 
/c^  by  degree  n,  and  also  cumulatively,  for  coefficient  set  B1  to  degree  250. 

The  same  information  is  also  listed  for  the  current  coefficient  set  A  (Rapp,  1981) 
till  degree  180. 

The  undulation,  and  anomaly  magnitude,  per  degree,  for  coefficient  set  B1 
have  been  plotted  in  Figures  6.1  and  6.2.  The  values  of  and  /c^  at  degree 
250  are  about  1  cm,  and  about  0.5  mgals  respectively. 


(rg)U|  UV%  ' 93j69q  jad  luauuaAOjduJi  a&ojuaojad 


Improvement  Due  to  Optimal  Estimation  of  Coefficients. 

Comparison  of  Both  Coefficients  Sets  A  and  B,  each  with  4  Optimally  Estimated  Coefficient  Sets 
from  Anomaly  Sets  A  and  B  Using  Error  Estimates  A  and  B  (Section  4). 
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anomaly  error  estimates  are  represented  by  error  estimate  A  in  Table  4.1.  Hence, 
a  slightly  conservative  estimate  of  current  improvement  in  coefficients  is  repre¬ 
sented  by  the  comparison  of  coefficient  B  (by  de-smoothing  factors  from  anomaly 
set  B)  with  the  optimally  estimated  coefficients  from  anomaly  set  B  using  error 
estimate  A.  This  improvement,  %An  in  (5.1),  is  seen  in  Table  5.7  to  be  about  8%  at 
degree  60,  about  11%  at  degree  120,  and  rises  to  about  33%  at  degree  180.  This  im¬ 
provement  is  also  shown  in  Figure  5.2.  The  improvement  by  degree  in  geo id  undu¬ 
lation  difference  in  (5.4),  and  anomaly  difference  /ac0  in  (5.5)  is  shown  in 
Figure  5.3.  The  improvements  have  been  shown  only  to  degree  180.  The  cumulative 
undulation  difference,  and  the  cumulative  anomaly  difference,  to  degree  180  was 
1.51  meters  and  2.70  mgal  respectively. 


cumulative  undulation  difference 


(5.6) 


cumulative  anomaly  difference 


rm 

=  \  l  ACn 


(5.7) 


Table  5.6  Limits  of  Four  10°x  10°  Blocks  where  l°x  1°  Anomalies  were  set  to  Zero 
in  Anomaly  Set  A. 


# 

Location 

(91 ) 

*s  (&i) 

m 

Remarks 

1 

Central  Siberia 

m 

1 

95 

105 

N.  Hemisphere 

2 

Central  Africa 

Wfimk 

1 

25 

35 

N.  Hemisphere 

3 

Southwest  Africa 

-10(100) 

-20(119) 

10 

20 

S.  Hemisphere 

4 

Central  Andes 
(S.  America) 

-20(120) 

-30(129) 

285 

295 

S.  Hemisphere 

By  comparing  the  first  two  sets  of  columns  in  Table  5.5,  we  find  noticeable 
changes  in  the  coefficients  when  200  l°x  1°  anomalies  are  set  to  zero  in  the  northern 
hemisphere  in  the  anomaly  set  A.  But  the  last  two  sets  of  columns  in  Table  5.5 
show  that  the  setting  of  200  l°x  1°  anomalies  to  zero  in  the  southern  hemisphere 
have  the  predominant  effect.  This  is  due  to  the  average  anomaly  error  estimate 
associated  with  the  latitude  band  in  which  the  anomaly  is  set  to  zero.  With  refer¬ 
ence  to  Tables  5.6  and  4.1,  the  value  of  o-j  associated  with  the  blocks  in  the 
northern  hemisphere  range  from  16  to  14  and  10  to  7  mgals,  while  ranges  from  7  to  8 
and  4  to  5  mgals  in  the  blocks  in  the  southern  hemisphere. 

These  tests  highlight  the  importance  of  getting  the  best  possible  estimate 
(instead  of  zero)  for  the  global  anomalies,  and  assigning  realistic  standard  errors 
to  these  estimates.  The  error  estimates  of  l°x  1°  anomalies  are,  of  course,  averaged 
over  each  latitude  band  before  these  are  utilized  in  the  optimal  estimation  procedure. 

5.4  Improvement  in  Coefficients  with  Optimal  Estimation 

We  now  list  in  Table  5.7  the  comparison  of  both  coefficient  sets  A  and  B  with 
optimally  estimated  coefficients  to  degree  180  from  both  anomaly  sets  A  and  B,  using 
in  each  case  anomaly  error  estimate  A  as  well  as  anomaly  error  estimate  B.  By 
listing  these  comparisons  side  by  side,  we  first  note  that  only  small  differences 
exist  in  comparison  with  coefficient  set  A,  when  contrasted  with  comparisons  made 
with  coefficient  set  B,  at  lower  degrees  £  12;  coefficient  sets  A  and  B  are,  of 
course,  exactly  the  same  for  degrees  >  36  (except  for  some  coefficients  to  degree  48). 
We  also  note  the  slightly  better  agreement  of  coefficients  estimated  from  anomaly  set 
B,  when  contrasted  with  similar  cases  of  coefficients  estimated  from  anomaly  set  A, 
except  at  lower  degrees  £  24.  The  large  disagreement  at  degrees  2  and  3  of  optimally 
estimated  coefficients  from  anomaly  set  B  was  commented  upon  in  Section  5.2.  The  dis¬ 
agreement  at  degree  3  becomes  worse  for  the  anomaly  error  estimate  A,  as  compared  to 
anomaly  error  estimate  B,  because  of  the  strong  equatorial  assymmetry  of  error  esti¬ 
mate  A  in  Table  4.1. 

Besides  other  points  already  discussed  in  Sections  5.1  to  5.3,  the  improvement, 
particularly  at  higher  degrees,  in  magnitude  of  coefficients  estimated  optimally 
instead  of  using  the  de-smoothing  factors  primarily  results  due  to  the  consideration 
of  realistic  anomaly  error  estimates  in  optimal  estimation,  while  the  anomaly  error 
estimates  are  not  considered  when  de-smoothing  factors  are  used.  The  current  realistic 


Variation  of  Optimally  Estimated  Coefficients  with  Variation  in  Global  Anomaly  Data. 
Comparison  of  Coefficient  Set  A  with  4  Coefficient  Sets  Estimated  from  4  Global  Anomaly  Sets 
Anomaly  Error  Estimate  A  (Table  4.1)  was  used  in  all  4  cases. 
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The  total  coefficient  percentage  error,  %Tn,  per  degree  has  been  plotted 
in  Figure  6.3  for  both  coefficient  sets  B1  and  B2.  Figure  6.3  has  been  repro¬ 
duced  from  Rapp  (1981,  p.  33,  Fig.  11)  which  showed  %Tn  for  the  current  coef¬ 
ficient  set  A.  %Tn  had  already  reached  100%  for  coefficient  set  A  at  degree 
120.  %Tn  is  better  for  sets  B1  and  B2  by  more  than  a  factor  of  two  as  compared 
to  set  A  for  all  degrees  greater  than  12. 


We  may  also  examine  the  error  estimates  in  undulation,  otn,  and  anomaly, 
acn,  per  degree  defined  as: 


=  R2  '  a2.„  =  R2  •  L_  (4nm  +  4J 


en 


m-0 


cn  =  y2(n-l)2  -«2£n  -  r2(n-l)2  l  (<&. m  +  4„J 


m=0 


ZTnm  "Srim 


(6.4) 

(6.5) 


and  also  cumulatively  as  ^a^zn  and  cn  respectively.  The  undulation  and 
anomaly  error  estimates,  per  degree,  and  cumulatively,  have  been  tabulated  in 
Table  6.4  for  coefficient  sets  Bl,  B2,  and  also  for  the  current  coefficient  set  A. 

We  note  that  error  estimates  for  sets  Bl  and  B2  are  better  than  a  factor  of  two 
from  the  corresponding  error  estimates  of  set  A  for  all  degrees  greater  than 
12.  But  the  error  estimates  for  sets  Bl  and  B2  are  much  larger  than  those  for 
set  A  at  degree  6  and  below,  though  as  noted  in  Table  6.3,  the  coefficient  per¬ 
centage  error  itself  for  sets  Bl  and  B2  is  quite  low  at  these  degrees.  This 
illustrates  the  well-known  fact  that  the  low  degree  potential  coefficients  can 
be  estimated  much  better  from  satellite  observations  as  compared  to  global 
anomaly  data.  The  error  estimates  of  low  degree  coefficients  in  set  A  have  been 
put  equal  to  the  error  estimates  of  '‘justed  SET1 ’  coefficients  (see  Section 
2.3),  while  the  error  estimates  of  sets  Bl  and  B2  are  based  on  the  anomaly  error 
estimates.  If  we  were  to  assume  that  the  error  estimates  for  set  Bl  for  degree 
£6  are  taken  from  satellite  observations,  i.e.  the  same  as  for  set  A,  then  the 
cumulative  undulation  error  for  set  Bl  to  degree  250  would  become  about  81  cm 
instead  of  109  cm.  In  a  similar  way,  the  cumulative  undulation  error  for  set 
B2  to  degree  250  would  be  reduced  to  about  56  cm  instead  of  67  cm.  The  cumulative 
undulation  errors  to  degree  180  may  then  be  compared  for  sets  Bl,  B2  and  A  as 
about  75,  48  and  146  cm  respectively.  The  error  estimates  in  Table  6.4,  and 
in  the  corresponding  Figures  6.4  and  6.5,  were  however  not  changed. 

The  cumulative  undulation  error,  ^°/£n ,  has  been  plotted  for  sets  A,  Bl 
and  B2  as  the  top  three  curves  in  Figure  6.4.  The  undulation  error,  a£n,  per 
degree  has  also  been  plotted  for  set  Bl  as  the  bottom  curve  in  Figure  6.4,  except 
that  the  portion  for  degrees  3,  4,  5  has  not  been  drawn  in  to  avoid  confusion. 

The  cumulative  anomaly  error,  '/5°2cn  ,  has  been  plotted  for  the  three  sets 
A,  Bl,  B2  in  Figure  6.5,  which  also  snows  the  anomaly  error,  acn,  per  degree  for 
set  Bl  as  the  bottom  curve. 


Table  6.4  Undulation  and  Anomaly  Error  Estimates  per  Degree  and  Cumulatively. 

Comparison  of  New  and  Current  Sets  of  Potential  Coefficients. 

Also,  Variation  Due  to  Anomaly  Error  Estimate  on  Optimally  Estimated  Coefficient  Accuracies. 
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Coefficient  set  A  is  the  current  set  of  coefficients  (Rapp,  1981)  to  degree  180.  (Also  see  Sec.  2.3).  o =  20  mgals, 
but  error  of  low  degree  coefficients  based  on  satellite  observations. 


6.3  Empirical  Anomaly  Degree  Variances 

It  is  of  interest  to  see  what  is  the  empirical  anomaly  degree  variance  cn, 
implied  by  the  optimally  estimated  coefficient  sets  B1  and  B2. 


= 


( n- 1 ) 2 


(6.6) 


Four  sets  of  anomaly  degree  variances  have  been  plotted  in  Figure  6.6.  At  the 
right  edge  of  the  plot,  the  curves  in  decreasing  order  of  magnitude  are:  (1)  the 
a-priori  cn  model  specified  in  (4.1)  plotted  with  asterisks;  (2)  Kaula's  model 
c|<  corresponding  to  =  10-5/n2  (Rapp,  1979,  p.  1,  (2)),  which  is  the  smooth 

curve  in  the  plot;  (3)  cn  for  coefficient  set  B2;  and  (4)  cn  for  coefficient 
set  B1  plotted  with  crosses.  As  the  optimal  estimator  in  (2.19)  in  least  squares 
collocation  minimizes  the  sum  of  total  variances  for  all  coefficients,  the  larger 
the  anomaly  error  estimate  the  smoother  would  be  the  optimally  estimated  coefficients. 
We  accordingly  find  that  cn  for  coefficient  set  B2  based  on  anomaly  estimate  B, 
aj  =  5  mgals,  have  consistently  larger  power  as  compared  to  coefficient  set  Bl, 
based  on  larger  anomaly  error  estimate  A,  with  ai  ranging  from  30  to  4  mgals 
(Table  4.1). 

We  also  note  that  the  power  in  coefficient  set  B2  falls  below  Kaula's  model 
for  degrees  higher  than  200,  while  the  a-priori  cn  model  has  too  much  power  for 
degrees  higher  than  150.  These  are  indicative  of  both  a  need  for  slight  downward 
adjustment  of  power  in  a-priori  model  in  (4.1),  and  also  that  more  reliable  esti¬ 
mates  for  higher  degrees  would  be  obtained  if  the  global  anomaly  data  was  available 
in  blocks  smaller  than  l°xl°. 

Because  of  greater  power  in  the  a-priori  model  for  cn,  the  error  estimates 
per  degree  are  pessimistic  at  higher  degrees,  as  larger  error  estimates  are  obtained 
from  (3.8)  with  larger  values  of  cn.  We  have  thus  pessimistic  undulation  and 
anomaly  error  estimates  per  degree,  azn  and  acn  in  Section  6.2,  for  higher  degrees. 

On  the  other  hand,  the  undulation  and  anomaly  magnitude  per  degree,  /j^  and  v'cfj’ 
in  Section  6.1,  decrease  rapidly  at  higher  degrees  for  coefficient  set  Bl  due  to 
large  anomaly  error  estimate  A.  This  explains  why  the  signal  to  noise  ratio 
appears  to  fall  below  1  for  coefficient  set  Bl  for  degrees  180  and  higher,  when 
we  compare  Tables  6.2  and  6.4.  Similar  results  occur  for  coefficient  set  B2, 
which  has  more  power  than  set  Bl,  for  degrees  200  and  higher.  This  artifact 
however  does  not  occur  in  Table  6.3  and  Figure  6.3,  where  the  coefficient  percen¬ 
tage  error  does  not  exceed  85%  even  at  degree  250.  This  is  due  to  the  percentage 
errors  being  computed  through  (6.2)  and  (6.3),  where  the  numerator  and  denominator 
both  depend  on  the  a-priori  cn  model. 

We  next  compare  £p  implied  by  coefficient  sets  Bl  and  B2,  with  cn  implied 
by  the  current  coefficient  set  (Rapp,  1981),  i.e.  coefficient  set  A  computed 
with  de-smoothing  factors  (2.7)  but  we  now  consider  the  latter  set  also  up  to 
degree  250,  including  the  sharp  discontinuity  at  degree  180.  cn  values  are  com¬ 
pared  for  coefficient  sets  B2  and  A  in  Figure  6.7;  and  for  coefficient  sets  Bl 
and  A  in  Figure  6.8.  The  curves  for  coefficient  set  A  are  plotted  with  asterisks 
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Figure  6.6  Anomaly  Degree  variances  c 

Implied  by:  (1)  A-priori  Model  (4.1) 
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in  both  figures,  while  the  curves  for  coefficient  sets  B2  and  B1  are  plotted 
with  crosses.  The  scale  of  cn  in  Figures  6.7  and  6.8  is  exaggerated  five  times 
as  compared  to  Figure  6.6. 

The  sharp  discontinuity  at  degree  180  for  coefficient  set  A  is  clearly  seen 
in  Figures  6.7  and  6.8.  The  smaller  discontinuity  at  degree  60  may  be  inferred 
in  Figure  6.7,  where  the  curve  with  asterisks  (coefficient  set  A)  has  more  power 
for  n<60,  and  lesser  power  for  n  >  60  as  compared  to  curve  with  crosses  (coef¬ 
ficient  set  B2).  Except  for  the  sharp  discontinuity  at  degree  180,  the  spectrums 
of  coefficient  sets  A  and  B2  are  fairly  close  in  Figure  6.7,  i.e.  when  we  do 
not  consider  any  latitudinal  variation  in  anomaly  error  estimate  B  for  the  optimal 
estimation  of  coefficient  set  B2.  However,  when  we  do  consider  the  currently 
realistic  latitudinal  variations  in  anomaly  error  estimate  A  for  the  optimal  esti¬ 
mation  of  coefficient  set  Bl,  the  spectrums  of  sets  B1  and  A  differ  substantially 
in  Figure  6.8  for  the  entire  range. 
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Figure  6.8  Comparison  of  Anomaly  Degree  Variances: 

Sets  A  and  B1 

Coefficient  Set  A  by  de-smoothing  factors  marked  by  asterisks. 

Coefficient  Set  B1  by  optimal  estimation  marked  by  crosses. 
Error  Estimate  A  :  a ^  =  30  to  4  mgals  (Table  4.1). 
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7.  SUMMARY  AND  CONCLUSIONS 


The  present  study  extends  the  algorithm  developed  by  Colombo  (1981),  and 
tested  with  5°x  5°  anomalies,  to  the  use  of  l°x  1°  anomalies  globally  for  the 
optimal  estimation  of  a  high  degree  potential  coefficients  field,  complete  to 
degree  and  order  250,  along  with  their  error  estimates.  The  algorithm  was  divided 
into  several  steps,  for  ease  of  computer  implementation  and  running  various  tests, 
and  also  needed  two  main  modifications.  These  steps  have  been  summarized  in 
Table  3.4. 

The  initial  first  step  is  to  assemble  the  integrals  of  associated  Legendre 
functions  in  1°  latitude  bands  by  an  efficient  algorithm  to  a  degree  Nmax,  which 
should  be  larger  than  the  highest  degree  NN  to  which  coefficients  will  be  esti¬ 
mated  and  smaller  than  (2ir/anomaly  block  size),  i.e.  360.  A  total  of 
(Nmax+1)(Nmax+2)/2  integrals,  which  comes  to  45,451  for  Nmax=300,  need  to  be 
computed  for  each  1°  latitude  band.  The  integrals  are  modified  by  factors,  de¬ 
pending  on  the  anomaly  degree  variances  and  the  anomaly  block  size,  for  later 
processing.  The  anomaly  degree  variance  model  in  (4.1)  was  used,  which  is 
the  current  global  model. 

The  generation  of  normals,  which  are  based  on  combinations  of  the  Fourier 
transforms  of  the  modified  integrals  in  the  first  step  for  different  pairs  of 
latitude  bands,  could  take  large  CPU  time  depending  on  degree  NN.  The  two  main 
modifications  to  the  algorithm,  i.e.  the  reordering  of  the  normals  and  the  re¬ 
ordering  of  the  right  hand  side  vectors,  became  necessary  due  to  very  large 
TAPE  10  and  array  size  requirements  for  the  case  of  l°xl°  anomalies.  These  two 
modifications,  which  are  described  in  Sections  3.31  and  3. 32, also  depend  on  degree 
NN.  As  none  of  the  steps  so  far  depend  on  anomaly  errors,  or  anomalies,  different 
tests  were  carried  out  as  described  in  Section  5  with  NN=12  or  60. 

The  solution  of  normals,  yielding  the  optimal  quadrature  weights,  and  the 
computation  of  error  estimates,  depends  on  the  anomaly  error  estimates.  An 
averaged  variance  is  used  in  each  latitude  band,  based  on  the  anomaly  error 
estimates  in  that  latitude  band,  as  discussed  in  Section  2.2  to  allow  the  inver¬ 
sion  of  covariance  matrix  of  64,800  l°x  1°  anomalies  to  be  equivalently  carried 
out  by  the  inversion  of  180  x  180  matrices  in  the  frequency  domain  to  degree  NN. 

The  use  of  averaged  anomaly  variance  in  each  latitude  band  allows  the  correct 
computation  of  only  the  sum  of  variances  for  the  two  coefficients  (Cnm,  Snm) 
for  any  particular  degree  and  order,  and  not  separate  variance  for  each  coefficient 
in  the  pair.  The  variance  for  each  coefficient  in  the  pair  is  then  arbitrarily 
made  equal . 


A  very  realistic  anomaly  error  estimate,  given  in  Table  4.1,  was  used 
for  solution  to  degree  NN-250.  Because  only  an  averaged  variance  is  required 
in  each  latitude  band,  the  anomaly  error  model  is  primarily  sensitive  to  latitude 
variation.  Accordingly,  another  solution  of  the  normals  was  done  to  degree 
NN=250  with  an  anomaly  error  estimate  without  any  latitude  variation.  The  two 
anomaly  error  estimates,  one  in  Table  4.1  with  average  anomaly  error  in  a  lati¬ 
tude  band  ranging  from  30  to  4  mgals,  and  the  other  with  a  uniform  value  of 
5  mgals  in  all  latitude  bands,  were  called  anomaly  error  estimates  A  and  B. 


The  computation  of  potential  coefficients  is  the  last  step,  where  the  global 
omaly  data  is  used  with  the  optimal  quadrature  weights,  which  are  dependent 
the  anomaly  error  estimates.  The  current  best  estimate  of  a  global  l°x  1° 
omaly  field,  anomaly  set  B  adjusted  to  a  weighted  set  of  satellite  derived 
tential  coefficients  ' SET1 1  complete  to  degree  36,  with  some  additional  coef- 
cients  to  degree  48,  was  described  in  Section  2.3.  Anomaly  set  B  was  used 
th  anomaly  error  estimates  A  and  B  to  compute  potential  coefficient  sets  B1 
d  B2  complete  to  degree  and  order  250.  The  statistics  on  these  two  coefficient 
ts,  and  the  comparison  with  the  current  set  of  potential  coefficients,  coef- 
cient  set  A,  complete  to  degree  and  order  180  (Rapp,  1981),  was  presented  in 
ction  6. 

When  the  next  updated  version  of  a  global  set  of  l°x  1°  anomalies  is  assembled, 
will  be  necessary  to  only  repeat  the  last  step  of  the  computation  of  an  up- 
ited  set  of  potential  coefficients.  The  latitude  variations  of  the  averaged 
lomaly  variances  in  anomaly  error  estimate  A  are  unlikely  to  change  in  the  near 
iture  till  the  flight  of  the  Geopotential  Research  Mission.  The  optimal  quadrature 
sights  corresponding  to  anomaly  error  estimate  A,  (and  those  corresponding  to 
lomaly  estimate  B  with  no  latitude  variation  of  anomaly  errors),  of  this  study 
lould  therefore  be  directly  usable  with  an  updated  version  of  global  l°x  1°  anomalies. 

.1  Conclusions  of  Present  Study 

-  The  most  reliable  estimate  should  be  used  for  the  global  gravity  anomaly 
sta.  When  no  reliable  l°x  1°  anomaly  estimates  a^e  available,  the  anomaly  implied 
/  a  high  degree  potential  coefficients  field  should  be  used  in  preference  to 
jtting  the  anomaly  estimate  as  zero.  The  lower  the  averaged  variance  of  anomalies 
i  a  1°  latitude  band,  the  larger  is  the  effect  of  change  in  anomaly  estimate  for 
vy  l°x  1°  anomaly  in  that  latitude  band. 

-  The  most  realistic  error  estimate  should  be  used  for  each  anomaly,  instead 
F  using  a  global  average.  Though  the  anomaly  variances  are  averaged  in  each 

3  latitudinal  band  before  their  introduction  in  optimal  estimation  of  coefficients, 
le  potential  coefficient  determination  is  quite  sensitive  to  any  resulting  latitude 
iriation  of  averaged  variances  in  each  1°  latitude  band. 

-  The  optimal  quadrature  weights,  Xnm,  are  computed  as  180  numbers  for  each 
secific  coefficient  pair,  (Cnm,  Snm),  of  given  degree  and  order.  The  180  numbers 
jke  into  account  the  latitudinal  variation  of  averaged  anomaly  variances.  The 
‘-smoothing  factors  are  based  on  the  Pel  1 inen/Meissl  smoothing  operator,  3n, 
lich  is  the  same  for  all  orders  for  a  given  degree.  Further,  the  de-smoothing 
jctors  do  not  take  account  of  any  anomaly  error  estimates,  or  more  specifically 
iy  latitudinal  variation  in  these  estimates. 

-  The  optimal  quadrature  weights  are  computed  based  on  the  minimization  of 
ie  propagated  anomaly  error  variances,  and  the  error  variance  due  to  sampling 

finite  number  of  mean  anomalies,  instead  of  an  infinite  number  of  point 
nomalies.  The  desmoothing  factors  are  based  only  on  the  size  of  anomaly  blocks. 

-  The  order  of  improvement  in  the  coefficients,  per  degree,  by  optimal 
stimation  over  the  current  estimates  using  de-smoothing  factors,  is  about  8% 
t  degree  60,  about  11%  at  degree  120,  and  rises  to  about  33%  at  degree  180. 


-  The  error  estimates  in  optimal  estimation  are  consistent  with  anomaly 
•ror  estimates  in  a  meaningful  way.  The  total  percentage  error,  %Tn,  per  degree 
les  not  exceed  100%  even  at  degree  250.  %Tn  had  exceeded  100%  for  current  set 

r  coefficients  at  degree  120. 

-  The  error  estimates  per  degree  in  optimal  estimation  are  better  than  a 
ictor  of  two  over  the  corresponding  estimates  in  the  current  set  of  coefficients 
>r  all  degrees  greater  than  12.  The  error  estimates  for  low  degrees  in  current 
>efficients  are  based  on  satellite  observations,  while  the  error  estimates  per 
:gree  in  optimal  estimation  are  based  only  on  anomaly  data. 

-  If  the  error  estimates  per  degree  in  optimal  estimation  procedure  are  set 
;  the  same  values  as  current  coefficient  set  only  for  degrees  £  6,  the  cumula- 
ive  undulation  error  to  degree  180  would  be  75  cm  as  compared  to  146  cm  for 
jrrent  set  of  coefficients. 

-  There  is  no  discontinuity  in  the  degree  variances  of  optimally  estimated 
efficients  at  degrees  60  and  180,  as  is  the  case  with  the  current  set  of  co- 
fficients,  which  were  limited  to  degree  180  because  of  concern  over  sharp  dis- 
jntinuity  at  that  degree. 

-  No  useful  purpose  may  be  served  by  expanding  l°x  1°  gravity  anomaly  data 
ito  coefficients  beyond  degree  250.  The  undulation  and  anomaly  magnitude  per 
agree  near  n=250  is  about  1  cm  and  0.5  mgals  respectively.  The  coefficient 
;rcentage  error  near  n=250  exceeds  80%. 

-  The  a-priori  anomaly  degree  variance  model  has  too  much  power  beyond  degree 
30.  This  leads  to  pessimistic  error  estimates  per  degree,  and  cumulatively,  at 
igher  degrees. 


.2  Recommendations  for  Further  Investigations 

-  The  optimal  estimation  algorithm  applied  to  l°x  1°  anomaly  data  needs  to 
e  extended  further  to  consider  the  case  when  global  anomaly  data  may  be  avail- 
ble  in  30'x30'  blocks.  Efforts  have  already  been  initiated  (private  communication 
rom  R.H.  Rapp)  to  assemble  a  global  data  of  259,200  30'x30'  anomalies.  This 

ill  lead  to  several  fold  increase  in  the  complexity  of  computer  implementation 
f  the  optimal  estimation  algorithm. 

-  The  generation  of  normals  requires  (Nmax+l)(Nmax+2)/2  integrals  of  associated 
egendre  functions  to  be  read  in  for  each  latitude  band  to  obtain  Fourier  transform 
or  the  data  in  a  pair  of  latitude  bands.  Even  after  exploiting  the  persymmetric 
tructure  of  the  normals  matrix,  about  136,000  TAPE  10 ' s  are  required  only  to 

ead  in  the  tiatefromall  required  pairs  of  latitude  bands  for  Nmax=300  for  1° 
atitude  bands.  It  needs  to  be  investigated  if  data  from  several  pairs  of  latitude 
ands  could  be  operated  upon  simultaneously  instead  of  one  pair  at  a  time.  Ef- 
iciency  of  putting  part  of  this  data  on  random  storage  also  needs  to  be  considered. 

-  The  CPU  time  for  the  solution  of  normals  increases  sharply  with  increase 
n  the  highest  degree  to  which  coefficients  are  estimated.  At  degree  250,  the 
olution  of  normals  in  sequence  of  right  hand  side  vectors  after  inverting 
atrices  of  only  180  x  180  for  the  case  of  l°x  1°  anomalies  already  took  about 

4  minutes  CPU  time  on  an  Amdahl  470  V/8  computer.  Faster  algorithms  need  to 
e  investigated  for  this  solution. 


-  Bose  et  al .  (1983)  have  suggested  some  strategies  for  including  data  at 
soles  but  the  computational  load  does  not  decrease  significantly.  They  have 
suggested  different  data  sampling  instead  of  equi-angular  blocks.  Further 
stigations  need  to  be  continued  to  apply  such  strategies  to  actual  computations 
high  degree  field. 


■ ■ -i 


-58- 


.  ■ i  *y  *  i v  ~  i  1  ',1.1,1  .'ii,  \  ,  V ■  >  \  i 


References 

,  S.C.,  G.E.  Thobe,  J.T.  Kouba  and  R.W.  Mortensen,  Optimal  Global  Gravity 
Representation,  Paper  presented  at  the  XVIII  General  Assembly  of  IUGG,  Proc. 
of  IAG  Symposia,  Vol .  1,  pp.  449-482,  Hamburg,  August  1983. 

nbo,  O.L.  ,  Optimal  Estimation  from  Data  Regularly  Sampled  on  a  Sphere  with 
Applications  in  Geodesy,  Department  of  Geodetic  Science  Report  No.  291, 

The  Ohio  State  University,  Columbus,  September  1979. 

nbo,  O.L.  ,  Numerical  Methods  for  Harmonic  Analysis  on  the  Sphere,  Department 
of  Geodetic  Science  and  Surveying,  Report  NO.  310,  The  Ohio  State  University, 
Columbus,  March  1981. 

tl ,  M. ,  On  the  Recursive  Computation  of  the  Integrals  of  the  Associated  Legendre 
Functions,  Manuscripta  Geodaetica,  Vol.  5,  pp.  181-199,  1980. 

son,  D.M.,  On  Solving  the  Stability  Problem  in  Computing  the  Integrals  of 
the  Legendre  Functions,  Department  of  Geodetic  Science  and  Surveying  Internal 
Report,  The  Ohio  State  University,  Columbus,  August  1983. 

h,  F.J.,  B.J.  Putney,  C.A.  Wagner  and  S.M.  Klosko,  Goddard  Earth  Models  for 
Oceanographic  Applications  (GEM  10B  and  10C),  Marine  Geodesy,  Vol.  5,  No.  2, 
pp.  145-187,  1981. 

tz,  H.,  Advanced  Physical  Geodesy,  Abacus  Press,  Kent,  U.K.,  1980a. 

tz,  H.,  Geodetic  Reference  System  1980,  Bulletin  Geodesique,  Vol.  54,  No.  3, 
pp.  395-405,  1980b. 

,  M.K.,  Recurrence  Relations  for  Integrals  of  Associated  Legendre  Functions, 
Bulletin  Geodesique,  Vol.  52,  pp.  177-190,  1978. 

i,  R.H.,  Analytical  and  Numerical  Differences  Between  Two  Methods  for  the 
Combination  of  Gravimetric  and  Satellite  Data,  Boll,  di  Geof.  Teo.  ed  Appl . , 
Vol.  XI,  N41-42,  pp.  108-118,  1979. 

i,  R.H.,  The  Relationship  Between  Mean  Anomaly  Block  Sizes  and  Spherical 
Harmonic  Representations,  Journal  of  Geophysical  Research,  Vol.  82,  No.  33, 
pp.  5360-5364,  November  1977. 

i,  R.H.,  Potential  Coefficient  and  Anomaly  Degree  Variance  Modeling  Revisited, 
Department  of  Geodetic  Science  and  Surveying  Report  No.  293,  The  Ohio  State 
University,  Columbus,  September  1979. 

i,  R.H.,  f he  Earth's  Gravity  Field  to  Degree  and  Order  180  Using  Seasat 
Altimeter  Data,  Terrestrial  Gravity  Data,  and  Other  Data,  Department  of 
Geodetic  Science  and  Surveying  Report  No.  322,  The  Ohio  State  University, 
Columbus,  December  1981. 

),  R.H.,  Aspects  of  Geoid  Definition  and  Determination,  Proc.  of  the  General 
Meeting  of  the  IAG,  pp.  411-433,  Tokyo,  May  1982. 


-59- 


Tscherning,  C.C.,  The  Role  of  High  Degree  Spherical  Harmonic  Expansions  in  Solving 
Geodetic  Problems,  Paper  presented  at  the  XVIII  General  Assembly  of  IUGG, 
Proc.  of  IAG  Symposia,  Vol.  1,  pp.  431-441,  Hamburg,  August  1983. 

Tscherning,  C.C.,  R.H.  Rapp  and  C.C  Goad,  A  Comparison  of  Methods  for  Computing 
Gravimetric  Quantities  from  High  Degree  Spherical  Harmonic  Expansions, 
Manuscripta  Geodaetica,  Vol.  8,  pp.  249-272,  1983. 


-60- 


END 

f- 

r 

I 

filmed 

•-85 


